Nothing
#' Function to calculate the test vector for an object fitted with \code{gamm4}
#'
#' @param mod an object fitted with \code{gamm4}
#' @param name character; name of the covariate for which inference should be
#' calculated
#' @param sigma2 variance to be used in the covariance definition. If \code{NULL},
#' the estimate \code{mod$gam$sig2} is used.
#' @param nrlocs number of locations at which p-values and intervals are to be computed
#' for non-linear terms. This directly corresponds to a sequence of \code{nrlocs} quantiles
#' of the given covariate values.
#' @param complete_effect logical; TRUE performs a (conservative) test whether
#' the whole spline has a significant influence after accounting for all other effects.
#' @details
#' Function provides the test vectors for every location of the given covariate
#'
#'
testvec_for_gamm4 <- function(mod, name, sigma2 = NULL, nrlocs=7,
complete_effect = FALSE)
{
if (grepl(",",name)) name <- unlist(strsplit(name, ","))
modmat <- model.matrix(mod$gam)
names <- attr(modmat, "dimnames")[[2]]
xnames <- attr(mod$gam$terms, "term.labels")
pnames <- c(name, setdiff(xnames, name))
# ind <- grep(name, names)
if("gamm" %in% class(mod)) {
sterms <- c(unlist(sapply(mod$gam$smooth,"[[","term")), names(mod$lme$coefficients$random))
} else {
sterms <- c(unlist(sapply(mod$gam$smooth,"[[","vn")), names(mod$mer@cnms))
}
SigmaInv <- mod$gam$Vp
if(is.null(sigma2)) sigma2 <- mod$gam$sig2
if(complete_effect){
# Adapted from iboost package (https://github.com/davidruegamer/iboost/)
# from the selectiveInference package
svdu_thresh <- function(x) {
svdx <- svd(x)
inds <- svdx$d > svdx$d[1] * sqrt(.Machine$double.eps)
return(svdx$u[, inds, drop = FALSE])
}
if (length(name) == 1) {
ind_effect <- which(grepl(paste0("s(",name), colnames(modmat), fixed = TRUE))
}
if (length(name) == 2) {
ind_effect <- which(grepl(paste0("te(",paste(name, collapse = ",")), colnames(modmat), fixed = TRUE))
}
Xj <- modmat[,ind_effect]
Xminusj <- svdu_thresh(modmat[,setdiff(1:ncol(modmat),ind_effect)])
vTs <- list(t(svdu_thresh(Xj - tcrossprod(Xminusj) %*% Xj)))
attr(vTs[[1]], "var") <- sigma2
}else{ # single point on a (non-)linear function
if(all(name %in% sterms))
{
if(length(nrlocs)>1) qs <- nrlocs else qs <-
quantile(mod$gam$model[[name]], seq(0,1,l=nrlocs+2)[-c(1,nrlocs+2)])
datap <- as.data.frame(cbind(qs, matrix(0, nrow = length(qs),
ncol = length(pnames) - 1)))
this_name <- paste0("s(",name,")")
}else{
datap <- as.data.frame(cbind(1, matrix(0, nrow = 1,
ncol = length(pnames) - 1)))
qs <- 1
this_name <- paste0("\\b",name,"\\b")
}
names(datap) <- pnames[apply(sapply(pnames, grepl, x = names), 2, any)]
lpm <- predict(mod$gam, newdata = datap, type = "lpmatrix")
# set columns to zero not associated with the variable
if(all(name %in% sterms)){
deactivate <- which(!grepl(this_name, colnames(lpm), fixed=TRUE))
}else{
deactivate <- which(!grepl(this_name, colnames(lpm)))
}
lpm[, deactivate] <- 0
if("gamm" %in% class(mod)){
hm_obj <- hatmatfun_gamm(mod)
vTs <- lapply(1:nrow(lpm), function(j){
k <- lpm[j,]%*%hm_obj[[2]]
attr(k, "var") <- as.numeric(lpm[j,]%*%hm_obj[[1]]%*%lpm[j,])*sigma2/mod$gam$sig2
attr(k, "loc") <- qs[j]
return(k)
})
}else{
warning("gamm4 inference currently not well tested.")
vTs <- lapply(1:nrow(lpm), function(j){
k <- lpm[j,]%*%SigmaInv%*%t(modmat)/mod$gam$sig2
attr(k, "var") <- as.numeric(lpm[j,]%*%SigmaInv%*%lpm[j,])*sigma2/mod$gam$sig2
attr(k, "loc") <- qs[j]
return(k)
})
}
}
return(vTs)
}
#' Function to calculate the test vector for an \code{merMod} object
#'
#' @param mod \code{merMod} object
#' @param which integer; if NULL, test vector is created for all
#' coefficients, otherwise for the which'th coefficient
#' @param marginal logical; whether to construct a test vector
#' from the marginal or condition perspective
#' @param VCOV covariance used for the construction of test vector
#' @param efficient logical; whether to compute the
#' test vector corresponding to the efficient beta estimator
#' in the marginal case
#' @param G true random effect covariance
#' @param sig2 true error variance
#'
#' @details
#' Function provides 2 different marginal test vectors
#' (efficient = TRUE / FALSE) and one conditional test vector
#'
#' Note that the covariance of residuals R is assumed to be diagonal.
#'
#'
#'
#'
testvec_for_mm <- function(
mod,
which = NULL,
marginal = TRUE,
VCOV = NULL,
G = NULL,
sig2 = NULL,
efficient = TRUE
)
{
if(class(mod)=="lm"){
X <- model.matrix(mod)
vTs <- solve(crossprod(X))%*%t(X)
if(is.null(which)) which <- 1:ncol(X)
return(lapply(which, function(j) t(vTs[j,])))
}
if(is.null(sig2)) sig2 <- sigma(mod)^2
Z <- getME(mod, "Z")
X <- getME(mod, "X")
n <- nrow(X)
Lambda <- getME(mod, "Lambda")
Lambdat <- getME(mod, "Lambdat")
if(!marginal){
if(length(mod@Gp)>3) stop("Please implement for more than 2 REs.")
C <- cbind(X,Z)
if(is.null(G) || dim(G)[2] != dim(getME(mod,"ST")[[1]])[2]){
A <- bdiag(list(matrix(0, ncol=ncol(X), nrow=ncol(X)), Lambda%*%Lambdat))
}else{
if(all(c(G)==0))
solveG <- G else solveG <- solve(G)
A <- bdiag(list(matrix(0, ncol=ncol(X), nrow=ncol(X)),
bdiag(list(solveG)[rep(1,ncol(Z)/ncol(G))])))
}
if(class(A)=="try-error") # G does not fit to the selected model
{
A <- bdiag(list(matrix(0, ncol=ncol(X), nrow=ncol(X)), Lambda%*%Lambdat))
warning("Could not create A, using the estimated covariance")
}
if(is.null(VCOV)) VCOV <- Matrix::crossprod(C)/sig2 + A
}
if(marginal){
L <- getME(mod, "L")
if(is.null(VCOV)){
V0inv <- diag(rep(1, n)) - crossprod(solve(L, system = "L") %*%
solve(L, Lambdat, system = "P") %*% t(Z))
# = solve(diag(n)+Z%*%Lambda%*%Lambdat%*%t(Z))
V0inv <- V0inv / sig2
}else{
V0inv <- solve(VCOV)
}
}
if(marginal){
if(efficient)
vTs <- solve(crossprod(X, V0inv%*%X))%*%t(X)%*%V0inv else
vTs <- solve(crossprod(X))%*%t(X)
}else{ # conditional
V0inv <- solve(VCOV)
vTs <- V0inv%*%Matrix::t(C)/sig2
}
if(is.null(which)) if(marginal) which <- 1:ncol(X) else which <- 1:ncol(C)
vTs <- lapply(which, function(j){
vt <- vTs[j,]
if(!marginal) attr(vt, "var") <- V0inv[j,j]
return(vt)
})
return(lapply(vTs,t))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.