1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
|
getSummary.mclogit <- function(obj,
alpha=.05,
rearrange=NULL,
...){
smry <- summary(obj)
N <- obj$N
coef <- smry$coefficients
varPar <- smry$varPar
lower.cf <- qnorm(p=alpha/2,mean=coef[,1],sd=coef[,2])
upper.cf <- qnorm(p=1-alpha/2,mean=coef[,1],sd=coef[,2])
coef <- cbind(coef,lower.cf,upper.cf)
colnames(coef) <- c("est","se","stat","p","lwr","upr")
if(length(rearrange)){
coef.grps <- lapply(rearrange,function(ii){
if(is.character(ii) && !all(ii %in% rownames(coef)))
stop("coefficient(s) ",dQuote(unname(ii[!(ii %in% rownames(coef))]))," do not exist")
structure(coef[ii,],
dimnames=list(names(ii),dimnames(coef)[[2]])
)
})
grp.titles <- names(rearrange)
coef.grps <- do.call(memisc::collect,coef.grps)
coef <- array(NA,dim=c(
dim(coef.grps)[1] + NROW(varPar),
dim(coef.grps)[2],
dim(coef.grps)[3]
))
coef[seq(dim(coef.grps)[1]),,] <- coef.grps
if(length(varPar))
coef[dim(coef.grps)[1]+seq(nrow(varPar)),,1] <- varPar
dimnames(coef) <- list(
c(dimnames(coef.grps)[[1]],rownames(varPar)),
dimnames(coef.grps)[[2]],
grp.titles
)
}
VarPar <- NULL
VarCov <- smry$VarCov
se_VarCov <- smry$se_VarCov
for(i in seq_along(VarCov)){
lv.i <- names(VarCov)[i]
vc.i <- VarCov[[i]]
vr.i <- diag(vc.i)
cv.i <- vc.i[lower.tri(vc.i)]
se_vc.i <- se_VarCov[[i]]
se_vr.i <- diag(se_vc.i)
se_cv.i <- se_vc.i[lower.tri(se_vc.i)]
nms.i <- rownames(vc.i)
nms.i <- gsub("(Intercept)","1",nms.i,fixed=TRUE)
vrnames.i <- paste0("Var(~",nms.i,"|",lv.i,")")
cvnames.i <- t(outer(nms.i,nms.i,FUN=paste,sep=":"))
cvnames.i <- cvnames.i[lower.tri(cvnames.i)]
if(length(cvnames.i))
cvnames.i <- paste0("Cov(~",cvnames.i,"|",lv.i,")")
vp.i <- matrix(NA,nrow=length(vr.i)+length(cv.i),ncol=6)
vp.i[,1] <- c(vr.i,cv.i)
vp.i[,2] <- c(se_vr.i,se_cv.i)
dimnames(vp.i) <- list(c(vrnames.i,cvnames.i),
c("est","se","stat","p","lwr","upr"))
VarPar <- c(VarPar,structure(list(vp.i),names=lv.i))
}
phi <- smry$phi
LR <- smry$null.deviance - smry$deviance
df <- obj$model.df
deviance <- deviance(obj)
if(df > 0){
p <- pchisq(LR,df,lower.tail=FALSE)
L0.pwr <- exp(-smry$null.deviance/N)
LM.pwr <- exp(-smry$deviance/N)
McFadden <- 1- smry$deviance/smry$null.deviance
Cox.Snell <- 1 - exp(-LR/N)
Nagelkerke <- Cox.Snell/(1-L0.pwr)
}
else {
LR <- NA
df <- NA
p <- NA
McFadden <- NA
Cox.Snell <- NA
Nagelkerke <- NA
}
ll <- obj$ll
AIC <- AIC(obj)
BIC <- AIC(obj,k=log(N))
sumstat <- c(
phi = phi,
LR = LR,
df = df,
#p = p,
logLik = ll,
deviance = deviance,
McFadden = McFadden,
Cox.Snell = Cox.Snell,
Nagelkerke = Nagelkerke,
AIC = AIC,
BIC = BIC,
N = N
)
ans <- list(coef= coef)
ans <- c(ans,VarPar)
if(length(smry$ngrps)){
G <-as.integer(smry$ngrps)
names(G) <- names(smry$ngrps)
names(G) <- paste("Groups by",names(G))
ans <- c(ans,list(Groups=G))
}
c(ans,
list(sumstat=sumstat,
call=obj$call,
contrasts = obj$contrasts,
xlevels = obj$xlevels))
}
getSummary.mmclogit <- getSummary.mclogit
|
