1################################################################################ 2#' 3#' @title Ajusted covariance matrix for linear mixed models according 4#' to Kenward and Roger 5#' @description Kenward and Roger (1997) describbe an improved small 6#' sample approximation to the covariance matrix estimate of the 7#' fixed parameters in a linear mixed model. 8#' @name kr-vcov 9#' 10################################################################################ 11## Implemented in Banff, august 2013; Søren Højsgaard 12 13#' @aliases vcovAdj vcovAdj.lmerMod vcovAdj_internal vcovAdj0 vcovAdj2 14#' vcovAdj.mer LMM_Sigma_G get_SigmaG get_SigmaG.lmerMod get_SigmaG.mer 15#' 16#' @param object An \code{lmer} model 17#' @param details If larger than 0 some timing details are printed. 18#' @return \item{phiA}{the estimated covariance matrix, this has attributed P, a 19#' list of matrices used in \code{KR_adjust} and the estimated matrix W of 20#' the variances of the covariance parameters of the random effetcs} 21#' 22#' \item{SigmaG}{list: Sigma: the covariance matrix of Y; G: the G matrices that 23#' sum up to Sigma; n.ggamma: the number (called M in the article) of G 24#' matrices) } 25#' 26#' @note If $N$ is the number of observations, then the \code{vcovAdj()} 27#' function involves inversion of an $N x N$ matrix, so the computations can 28#' be relatively slow. 29#' @author Ulrich Halekoh \email{uhalekoh@@health.sdu.dk}, Søren Højsgaard 30#' \email{sorenh@@math.aau.dk} 31#' @seealso \code{\link{getKR}}, \code{\link{KRmodcomp}}, \code{\link{lmer}}, 32#' \code{\link{PBmodcomp}}, \code{\link{vcovAdj}} 33#' @references Ulrich Halekoh, Søren Højsgaard (2014)., A Kenward-Roger 34#' Approximation and Parametric Bootstrap Methods for Tests in Linear Mixed 35#' Models - The R Package pbkrtest., Journal of Statistical Software, 36#' 58(10), 1-30., \url{https://www.jstatsoft.org/v59/i09/} 37#' 38#' Kenward, M. G. and Roger, J. H. (1997), \emph{Small Sample Inference for 39#' Fixed Effects from Restricted Maximum Likelihood}, Biometrics 53: 983-997. 40#' 41#' @keywords inference models 42#' @examples 43#' 44#' fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) 45#' class(fm1) 46#' 47#' ## Here the adjusted and unadjusted covariance matrices are identical, 48#' ## but that is not generally the case: 49#' 50#' v1 <- vcov(fm1) 51#' v2 <- vcovAdj(fm1, details=0) 52#' v2 / v1 53#' 54#' ## For comparison, an alternative estimate of the variance-covariance 55#' ## matrix is based on parametric bootstrap (and this is easily 56#' ## parallelized): 57#' 58#' \dontrun{ 59#' nsim <- 100 60#' sim <- simulate(fm.ml, nsim) 61#' B <- lapply(sim, function(newy) try(fixef(refit(fm.ml, newresp=newy)))) 62#' B <- do.call(rbind, B) 63#' v3 <- cov.wt(B)$cov 64#' v2/v1 65#' v3/v1 66#' } 67#' 68#' 69#' 70#' @export vcovAdj 71#' 72#' @rdname kr-vcov 73vcovAdj <- function(object, details=0){ 74 UseMethod("vcovAdj") 75} 76 77#' @method vcovAdj lmerMod 78#' @rdname kr-vcov 79#' @export vcovAdj.lmerMod 80vcovAdj.lmerMod <- function(object, details=0){ 81 if (!(getME(object, "is_REML"))) { 82 object <- update(object, . ~ ., REML = TRUE) 83 } 84 Phi <- vcov(object) 85 SigmaG <- get_SigmaG(object, details) 86 X <- getME(object, "X") 87 vcovAdj_internal(Phi, SigmaG, X, details=details) 88} 89 90## Needed to avoid emmeans choking. 91 92#' @export 93vcovAdj.lmerMod <- vcovAdj.lmerMod 94 95 96 97## Dette er en kopi af '2015' udgaven 98vcovAdj_internal <- function(Phi, SigmaG, X, details=0){ 99 100 details=0 101 DB <- details > 0 ## debugging only 102 t0 <- proc.time() 103 104 if (DB){ 105 cat("vcovAdj16_internal\n") 106 cat(sprintf("dim(X) : %s\n", toString(dim(X)))) 107 print(class(X)) 108 cat(sprintf("dim(Sigma) : %s\n", toString(dim(SigmaG$Sigma)))) 109 print(class(SigmaG$Sigma)) 110 } 111 112 113 ##SigmaInv <- chol2inv( chol( forceSymmetric(SigmaG$Sigma) ) ) 114 SigmaInv <- chol2inv( chol( forceSymmetric(as(SigmaG$Sigma, "matrix")))) 115 ##SigmaInv <- as(SigmaInv, "dpoMatrix") 116 117 if(DB){ 118 cat(sprintf("Finding SigmaInv: %10.5f\n", (proc.time()-t0)[1] )); 119 t0 <- proc.time() 120 } 121 122 #mat <<- list(SigmaG=SigmaG, SigmaInv=SigmaInv, X=X) 123 124 t0 <- proc.time() 125 ## Finding, TT, HH, 00 126 n.ggamma <- SigmaG$n.ggamma 127 TT <- SigmaInv %*% X 128 HH <- OO <- vector("list", n.ggamma) 129 for (ii in 1:n.ggamma) { 130 #.tmp <- SigmaG$G[[ii]] %*% SigmaInv 131 #HH[[ ii ]] <- .tmp 132 #OO[[ ii ]] <- .tmp %*% X 133 HH[[ ii ]] <- SigmaG$G[[ii]] %*% SigmaInv 134 OO[[ ii ]] <- HH[[ ii ]] %*% X 135 } 136 if(DB){cat(sprintf("Finding TT, HH, OO %10.5f\n", (proc.time()-t0)[1] )); t0 <- proc.time()} 137 138 ## Finding PP, QQ 139 PP <- QQ <- NULL 140 for (rr in 1:n.ggamma) { 141 OrTrans <- t( OO[[ rr ]] ) 142 PP <- c(PP, list(forceSymmetric( -1 * OrTrans %*% TT))) 143 for (ss in rr:n.ggamma) { 144 QQ <- c(QQ, list(OrTrans %*% SigmaInv %*% OO[[ss]] )) 145 }} 146 if(DB){cat(sprintf("Finding PP,QQ: %10.5f\n", (proc.time()-t0)[1] )); t0 <- proc.time()} 147 148 149 ##stat15 <<- list(HH=HH, OO=OO, PP=PP, Phi=Phi, QQ=QQ) 150 151 Ktrace <- matrix( NA, nrow=n.ggamma, ncol=n.ggamma ) 152 for (rr in 1:n.ggamma) { 153 HrTrans <- t( HH[[rr]] ) 154 for (ss in rr:n.ggamma){ 155 Ktrace[rr,ss] <- Ktrace[ss,rr]<- sum( HrTrans * HH[[ss]] ) 156 }} 157 if(DB){cat(sprintf("Finding Ktrace: %10.5f\n", (proc.time()-t0)[1] )); t0 <- proc.time()} 158 159 ## Finding information matrix 160 IE2 <- matrix( NA, nrow=n.ggamma, ncol=n.ggamma ) 161 for (ii in 1:n.ggamma) { 162 Phi.P.ii <- Phi %*% PP[[ii]] 163 for (jj in c(ii:n.ggamma)) { 164 www <- .indexSymmat2vec( ii, jj, n.ggamma ) 165 IE2[ii,jj]<- IE2[jj,ii] <- Ktrace[ii,jj] - 166 2 * sum(Phi * QQ[[ www ]]) + sum( Phi.P.ii * ( PP[[jj]] %*% Phi)) 167 }} 168 if(DB){cat(sprintf("Finding IE2: %10.5f\n", (proc.time()-t0)[1] )); t0 <- proc.time()} 169 170 eigenIE2 <- eigen(IE2, only.values=TRUE)$values 171 condi <- min(abs(eigenIE2)) 172 173 WW <- if (condi > 1e-10) forceSymmetric(2 * solve(IE2)) else forceSymmetric(2 * ginv(IE2)) 174 175 ## print("vcovAdj") 176 UU <- matrix(0, nrow=ncol(X), ncol=ncol(X)) 177 ## print(UU) 178 for (ii in 1:(n.ggamma-1)) { 179 for (jj in c((ii + 1):n.ggamma)) { 180 www <- .indexSymmat2vec( ii, jj, n.ggamma ) 181 UU <- UU + WW[ii,jj] * (QQ[[ www ]] - PP[[ii]] %*% Phi %*% PP[[jj]]) 182 }} 183 ## print(UU) 184 185 UU <- UU + t(UU) 186 ## UU <<- UU 187 for (ii in 1:n.ggamma) { 188 www <- .indexSymmat2vec( ii, ii, n.ggamma ) 189 UU<- UU + WW[ii, ii] * (QQ[[ www ]] - PP[[ii]] %*% Phi %*% PP[[ii]]) 190 } 191 ## print(UU) 192 GGAMMA <- Phi %*% UU %*% Phi 193 PhiA <- Phi + 2 * GGAMMA 194 attr(PhiA, "P") <-PP 195 attr(PhiA, "W") <-WW 196 attr(PhiA, "condi") <- condi 197 PhiA 198} 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241## #' @method vcovAdj mer 242## #' @rdname kr-vcov 243## #' @export 244## vcovAdj.mer <- vcovAdj.lmerMod 245 246 247 248## .vcovAdj_internal <- function(Phi, SigmaG, X, details=0){ 249 250## ##cat("vcovAdj_internal\n") 251## ##SG<<-SigmaG 252## DB <- details > 0 ## debugging only 253 254## #print("HHHHHHHHHHHHHHH") 255## #print(system.time({chol( forceSymmetric(SigmaG$Sigma) )})) 256## #print(system.time({chol2inv( chol( forceSymmetric(SigmaG$Sigma) ) )})) 257 258## ## print("HHHHHHHHHHHHHHH") 259## ## Sig <- forceSymmetric( SigmaG$Sigma ) 260## ## print("HHHHHHHHHHHHHHH") 261## ## print(system.time({Sig.chol <- chol( Sig )})) 262## ## print(system.time({chol2inv( Sig.chol )})) 263 264## t0 <- proc.time() 265## ## print("HHHHHHHHHHHHHHH") 266## SigmaInv <- chol2inv( chol( forceSymmetric(SigmaG$Sigma) ) ) 267## ## print("DONE --- HHHHHHHHHHHHHHH") 268 269## if(DB){ 270## cat(sprintf("Finding SigmaInv: %10.5f\n", (proc.time()-t0)[1] )); 271## t0 <- proc.time() 272## } 273## ##print("iiiiiiiiiiiii") 274 275## t0 <- proc.time() 276## ## Finding, TT, HH, 00 277## n.ggamma <- SigmaG$n.ggamma 278## TT <- SigmaInv %*% X 279## HH <- OO <- vector("list", n.ggamma) 280## for (ii in 1:n.ggamma) { 281## .tmp <- SigmaG$G[[ii]] %*% SigmaInv 282## HH[[ ii ]] <- .tmp 283## OO[[ ii ]] <- .tmp %*% X 284## } 285## if(DB){cat(sprintf("Finding TT,HH,OO %10.5f\n", (proc.time()-t0)[1] )); t0 <- proc.time()} 286## ## if(DB){ 287## ## cat("HH:\n"); print(HH); HH <<- HH 288## ## cat("OO:\n"); print(OO); OO <<- OO 289## ## } 290 291## ## Finding PP, QQ 292## PP <- QQ <- NULL 293## for (rr in 1:n.ggamma) { 294## OrTrans <- t( OO[[ rr ]] ) 295## PP <- c(PP, list(forceSymmetric( -1 * OrTrans %*% TT))) 296## for (ss in rr:n.ggamma) { 297## QQ <- c(QQ,list(OrTrans %*% SigmaInv %*% OO[[ss]] )) 298## }} 299## if(DB){cat(sprintf("Finding PP,QQ: %10.5f\n", (proc.time()-t0)[1] )); t0 <- proc.time()} 300## ## if(DB){ 301## ## cat("PP:\n"); print(PP); PP2 <<- PP 302## ## cat("QP:\n"); print(QQ); QQ2 <<- QQ 303## ## } 304 305## Ktrace <- matrix( NA, nrow=n.ggamma, ncol=n.ggamma ) 306## for (rr in 1:n.ggamma) { 307## HrTrans <- t( HH[[rr]] ) 308## for (ss in rr:n.ggamma){ 309## Ktrace[rr,ss] <- Ktrace[ss,rr]<- sum( HrTrans * HH[[ss]] ) 310## }} 311## if(DB){cat(sprintf("Finding Ktrace: %10.5f\n", (proc.time()-t0)[1] )); t0 <- proc.time()} 312 313## ## Finding information matrix 314## IE2 <- matrix( NA, nrow=n.ggamma, ncol=n.ggamma ) 315## for (ii in 1:n.ggamma) { 316## Phi.P.ii <- Phi %*% PP[[ii]] 317## for (jj in c(ii:n.ggamma)) { 318## www <- .indexSymmat2vec( ii, jj, n.ggamma ) 319## IE2[ii,jj]<- IE2[jj,ii] <- Ktrace[ii,jj] - 320## 2 * sum(Phi*QQ[[ www ]]) + sum( Phi.P.ii * ( PP[[jj]] %*% Phi)) 321## }} 322## if(DB){cat(sprintf("Finding IE2: %10.5f\n", (proc.time()-t0)[1] )); t0 <- proc.time()} 323 324## eigenIE2 <- eigen(IE2,only.values=TRUE)$values 325## condi <- min(abs(eigenIE2)) 326 327## WW <- if(condi>1e-10) forceSymmetric(2* solve(IE2)) else forceSymmetric(2* ginv(IE2)) 328 329## ## print("vcovAdj") 330## UU <- matrix(0, nrow=ncol(X), ncol=ncol(X)) 331## ## print(UU) 332## for (ii in 1:(n.ggamma-1)) { 333## for (jj in c((ii+1):n.ggamma)) { 334## www <- .indexSymmat2vec( ii, jj, n.ggamma ) 335## UU <- UU + WW[ii,jj] * (QQ[[ www ]] - PP[[ii]] %*% Phi %*% PP[[jj]]) 336## }} 337## ## print(UU) 338 339## UU <- UU + t(UU) 340## ## UU <<- UU 341## for (ii in 1:n.ggamma) { 342## www <- .indexSymmat2vec( ii, ii, n.ggamma ) 343## UU<- UU + WW[ii,ii] * (QQ[[ www ]] - PP[[ii]] %*% Phi %*% PP[[ii]]) 344## } 345## ## print(UU) 346## GGAMMA <- Phi %*% UU %*% Phi 347## PhiA <- Phi + 2 * GGAMMA 348## attr(PhiA, "P") <-PP 349## attr(PhiA, "W") <-WW 350## attr(PhiA, "condi") <- condi 351## PhiA 352## } 353