# R/predictFunctions.R In LMMsolver: Linear Mixed Model Solver

#### Documented in calcStandardErrorspredictTest

#' Standard errors for predictions
#'
#' Calculates the standard errors for predictions \eqn{D \hat{u}},
#' see Welham et al. 2004 and Gilmour et al. 2004 for details.
#'
#' The prediction error variance is given by \eqn{D C^{-1} D'},
#' where \eqn{C} is the mixed model coefficient matrix, and \eqn{D} defines
#' linear combinations of fixed and random effects.
#' The standard errors are given by the the square root of
#' the diagonal. To calculate the standard errors in an efficient way we use that
#'
#' \deqn{\frac{\partial log|C + \xi_i d_i d_i'|}{\partial \xi_i} |_{\xi_i=0}
#'  = trace(C^{-1} d_i d_i') =
#' trace(d_i' C^{-1} d_i) = d_i' C^{-1} d_i, }
#' where \eqn{d_i} is row \eqn{i} of matrix \eqn{D}. The values of
#' \eqn{d_i' C^{-1} d_i} can be calculated more efficient, avoiding the
#' calculation of the inverse of \eqn{C}, by using Automated Differentiation
#' of the Choleksy algorithm, see section 2.3 in Smith (1995) for details.
#'
#' @param C a symmetric matrix of class spam
#' @param D a matrix of class spam
#'
#' @return a vector with standard errors for predictions \eqn{D \hat{u}}.
#'
#' @references
#' Welham, S., Cullis, B., Gogel, B., Gilmour, A., & Thompson, R. (2004).
#' Prediction in linear mixed models.
#' Australian & New Zealand Journal of Statistics, 46(3), 325-347.
#'
#' @references
#' Smith, S. P. (1995). Differentiation of the Cholesky algorithm.
#' Journal of Computational and Graphical Statistics, 4(2), 134-147.
#'
#' @references
#' Gilmour, A., Cullis, B., Welham, S., Gogel, B., & Thompson, R. (2004).
#' An efficient computing strategy for prediction in mixed linear models.
#' Computational statistics & data analysis, 44(4), 571-586.
#'
#' @keywords internal
calcStandardErrors <- function(C,
D) {
## !!! NOT CHANGE THE LINE OF CODE BELOW !!!
## It adds extra zeros ("fill-ins") to matrix C, needed
## to calculate the Partial Derivatives of Cholesky, not equal to zero.
C = C + 0 * spam::crossprod.spam(D)
cholC <- chol(C)
p <- cholC@pivot
tD <- t(D)
tDp <- tD[p, ]
x <- diagXCinvXt(cholC, tDp)
se <- sqrt(x)
return(se)
}

auxFun <- function(object, i, nr, nc, Names=NULL)
{
# MB: More efficient to use t(Dg), working row-wise:
Dg <- spam::spam(x = 0, nrow = nr, ncol = nc)
ndx <- object$ndxCoefficients[[i]] len_ndx <- length(ndx) rm <- which(ndx==0) if (length(rm)>0) { ndx <- ndx[-which(ndx==0)] } if (!is.null(Names)) { for (j in 1:length(ndx)) { Dg[,ndx[j]] <- 1*(Names==names(ndx)[j]) } } else { # averaging Dg[,ndx] <- 1/len_ndx } Dg } #' Test function for predict, for the moment internal #' #' @keywords internal predictTest <- function(object, classify) { # count what kind of interactions, if any: cat("Predict function: \n") s <- classify s2 <- gsub(":","", s) numOcc <- nchar(s) - nchar(s2) numOcc if (numOcc == 1) { Interaction <- TRUE } else { Interaction <- FALSE } if (Interaction) { # names of marginals marg_names <- unlist(strsplit(classify, split=":")) } else { marg_names <- NULL } term.labels <- c(object$term.labels.f, object$term.labels.r) isFixed <- term.labels %in% object$term.labels.f
dim <- object$dim e <- cumsum(dim) s <- e - dim + 1 df.terms <- data.frame(label = term.labels, s = s, e = e, fixed = isFixed) ndxPred <- object$ndxCoefficients[[classify]]

namesI <- names(ndxPred)

nr <- length(ndxPred)
nc <- ncol(object$C) ndx_intercept <- which(df.terms$label == '(Intercept)')
ndx_classify <- which(df.terms$label == classify) ndx_marginals <- which(df.terms$label %in% marg_names)
ndx_averaging <- which(df.terms$fixed) ndx_averaging <- setdiff(ndx_averaging, ndx_intercept) ndx_averaging <- setdiff(ndx_averaging, ndx_classify) ndx_averaging <- setdiff(ndx_averaging, ndx_marginals) ndx_averaging cat(" Intercept: ", df.terms$label[ndx_intercept],"\n")
cat("  Classify:  ", df.terms$label[ndx_classify], "\n") cat(" Marginals: ", df.terms$label[ndx_marginals], "\n")
cat("  Averaging: ", df.terms$label[ndx_averaging], "\n") # intercept: Dg <- spam::spam(x = 0, nrow = nr, ncol = nc) Dg[, 1] <- 1 # ndx for classify: Dg <- Dg + auxFun(object, ndx_classify, nr, nc, namesI) for (i in ndx_marginals) { label <- df.terms$label[i]          # label of current term
k <- which(label==marg_names)
namesLab <- sapply(namesI, FUN = function(x) { strsplit(x,split=":")[[1]][k] })
Dg <- Dg + auxFun(object, i, nr, nc, namesLab)
}
for (i in ndx_averaging) {
Dg <- Dg + auxFun(object, i, nr, nc, NULL)
}
Dg <- spam::cleanup(Dg)

ypred <- as.vector(Dg %*% object$coefMME) ypredse <- calcStandardErrors(object$C, Dg)

if (Interaction) {
marg1_labels <- sapply(namesI, FUN = function(x) { strsplit(x,split=":")[[1]][1]})
marg2_labels <- sapply(namesI, FUN = function(x) { strsplit(x,split=":")[[1]][2]})

marg1 <- sapply(marg1_labels, FUN=function(x) { strsplit(x,"_")[[1]][2]})
marg2 <- sapply(marg2_labels, FUN=function(x) { strsplit(x,"_")[[1]][2]})
pred <- data.frame(marg1, marg2, ypred, ypredse)
colnames(pred) <- c(marg_names, "prediction", "se")
# order on first marginal
ord <- order(pred[[marg_names[1]]])
pred <- pred[ord, ]
} else {
term <- sapply(namesI, FUN=function(x) { strsplit(x,"_")[[1]][2]})
pred <- data.frame(term, ypred, ypredse)
colnames(pred) <- c(classify, "prediction", "se")
}

rownames(pred) <- NULL
pred
}


## Try the LMMsolver package in your browser

Any scripts or data that you put into this service are public.

LMMsolver documentation built on June 22, 2024, 12:16 p.m.