1getSummary.clm <- function(obj,
2                           alpha=.05,
3                           ...){
4
5  smry <- summary(obj)
6  N <- if(length(weights(obj))) sum(weights(obj))
7  else smry$nobs
8
9  prmtable <- coef(smry)
10
11  ## Move threshold parameters to end
12  coefnames <- names(smry$beta)
13  parnames <- rownames(prmtable)
14  threshnames <- setdiff(parnames,coefnames)
15
16  lower <- qnorm(p=alpha/2,mean=prmtable[,1],sd=prmtable[,2])
17  upper <- qnorm(p=1-alpha/2,mean=prmtable[,1],sd=prmtable[,2])
18
19  prmtable <- cbind(prmtable,lower,upper)
20
21  colnames(prmtable) <- c("est","se","stat","p","lwr","upr")
22
23  dn <- list(
24    rownames(prmtable),
25    c("est","se","stat","p","lwr","upr"),
26    names(obj$model)[1]
27  )
28  dim(prmtable) <- c(dim(prmtable)[1],dim(prmtable)[2],1)
29  dimnames(prmtable) <- dn
30
31  coef <- prmtable[coefnames,,,drop=FALSE]
32  thresh <- prmtable[threshnames,,,drop=FALSE]
33
34  ans <- list(coef = coef,
35              Thresholds = thresh)
36
37  null.model <- update(obj, .~1)
38
39  ll <- logLik(obj)
40  ll0 <- logLik(null.model)
41
42  LR <- 2*(ll-ll0)
43  df <- null.model$df.residual - smry$df.residual
44
45  dev <- -2*ll
46
47  if(df > 0){
48    p <- pchisq(LR,df,lower.tail=FALSE)
49    L0.pwr <- exp(2*ll0/N)
50
51    Aldrich.Nelson <- LR/(LR+N)
52    McFadden <- 1 - dev/(-2*ll0)
53    Cox.Snell <- 1 - exp(-LR/N)
54    Nagelkerke <- Cox.Snell/(1-L0.pwr)
55  }
56  else {
57    LR <- NA
58    df <- NA
59    p <- NA
60    Aldrich.Nelson <- NA
61    McFadden <- NA
62    Cox.Snell <- NA
63    Nagelkerke <- NA
64  }
65
66  AIC <- AIC(obj)
67  BIC <- AIC(obj,k=log(N))
68  sumstat <- c(
69    LR             = LR,
70    df         = df,
71    p             = p,
72    logLik        = ll,
73    deviance      = dev,
74    Aldrich.Nelson = Aldrich.Nelson,
75    McFadden      = McFadden,
76    Cox.Snell       = Cox.Snell,
77    Nagelkerke    = Nagelkerke,
78    AIC           = AIC,
79    BIC           = BIC,
80    N             = N
81  )
82
83  ans <- c(ans,
84           list(sumstat=sumstat,
85                contrasts=obj$contrasts,
86                xlevels=smry$xlevels,
87                call=obj$call))
88  return(ans)
89}
90
91getSummary.clmm <- function(obj,
92                           alpha=.05,
93                           ...){
94
95  smry <- summary(obj)
96  N <- if(length(weights(obj))) sum(weights(obj))
97       else nrow(obj$model)
98
99  prmtable <- coef(smry)
100
101  ## Move threshold parameters to end
102  coefnames <- names(smry$beta)
103  parnames <- rownames(prmtable)
104  threshnames <- setdiff(parnames,coefnames)
105
106  lower <- qnorm(p=alpha/2,mean=prmtable[,1],sd=prmtable[,2])
107  upper <- qnorm(p=1-alpha/2,mean=prmtable[,1],sd=prmtable[,2])
108
109  prmtable <- cbind(prmtable,lower,upper)
110
111  colnames(prmtable) <- c("est","se","stat","p","lwr","upr")
112
113  dn <- list(
114    rownames(prmtable),
115    c("est","se","stat","p","lwr","upr"),
116    names(obj$model)[1]
117  )
118  dim(prmtable) <- c(dim(prmtable)[1],dim(prmtable)[2],1)
119  dimnames(prmtable) <- dn
120
121  coef <- prmtable[coefnames,,,drop=FALSE]
122  thresh <- prmtable[threshnames,,,drop=FALSE]
123
124  ans <- list(coef = coef,
125              Thresholds = thresh)
126
127  varcor <- lapply(obj$ST,tcrossprod)
128
129  VarPar <- NULL
130
131  for(i in seq_along(varcor)){
132    vc.i <- varcor[[i]]
133    lv.i <- names(varcor)[i]
134    vr.i <- diag(vc.i)
135    cv.i <- vc.i[lower.tri(vc.i)]
136    nms.i <- rownames(vc.i)
137    nms.i <- gsub("(Intercept)","1",nms.i,fixed=TRUE)
138    vrnames.i <- paste0("Var(~",nms.i,"|",lv.i,")")
139    cvnames.i <- t(outer(nms.i,nms.i,FUN=paste,sep=":"))
140    cvnames.i <- cvnames.i[lower.tri(cvnames.i)]
141    if(length(cvnames.i))
142      cvnames.i <- paste0("Cov(~",cvnames.i,"|",lv.i,")")
143    vp.i <- matrix(NA,nrow=length(vr.i)+length(cv.i),ncol=6)
144    vp.i[,1] <- c(vr.i,cv.i)
145    dim(vp.i) <- c(dim(vp.i),1)
146    dimnames(vp.i) <- list(c(vrnames.i,cvnames.i),
147                           c("est","se","stat","p","lwr","upr"),
148                           names(obj$model)[1])
149    VarPar <- c(VarPar,structure(
150                           list(vp.i),
151                           names=lv.i))
152 }
153  ans <- c(ans,list(Variances=VarPar))
154
155  # null.model <- update(obj, .~1)
156
157  ll <- logLik(obj)
158  # ll0 <- logLik(null.model)
159
160  # LR <- 2*(ll-ll0)
161  # df <- null.model$df.residual - smry$df.residual
162
163  dev <- -2*ll
164
165  # if(df > 0){
166  #   p <- pchisq(LR,df,lower.tail=FALSE)
167  #   L0.pwr <- exp(2*ll0/N)
168  #   Aldrich.Nelson <- LR/(LR+N)
169  #   McFadden <- 1 - dev/(-2*ll0)
170  #   Cox.Snell <- 1 - exp(-LR/N)
171  #   Nagelkerke <- Cox.Snell/(1-L0.pwr)
172  # }
173  # else {
174    df <- NA
175    LR <- NA
176    p <- NA
177    Aldrich.Nelson <- NA
178    McFadden <- NA
179    Cox.Snell <- NA
180    Nagelkerke <- NA
181  # }
182
183  AIC <- AIC(obj)
184  # BIC <- AIC(obj,k=log(N))
185  sumstat <- c(
186    # LR            = LR,
187    # df            = df,
188    # p             = p,
189    logLik        = ll,
190    deviance      = dev,
191    # Aldrich.Nelson = Aldrich.Nelson,
192    # McFadden      = McFadden,
193    # Cox.Snell       = Cox.Snell,
194    # Nagelkerke    = Nagelkerke,
195    AIC           = AIC,
196  #   BIC           = BIC,
197    N             = N
198  )
199
200  ans <- c(ans,
201           list(sumstat=sumstat,
202                contrasts=obj$contrasts,
203                xlevels=smry$xlevels,
204                call=obj$call))
205  return(ans)
206}
207
208