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