varComp: Fitting variance component models

Description

varComp and varComp.fit fit linear mixed-effect models where the marginal variance-covariance matrix is linear in known positive semidefinite matrices. varComp uses the usual formula interface, whereas varComp.fit is the underlying working horse. doFit.varComp performs model fitting if the object has been previously created by setting doFit = FALSE when calling varComp.

Usage

varComp(fixed, data, random, varcov, weights, subset, 
family = stats::gaussian('identity'), na.action, offset, 
control = varComp.control(...), doFit = TRUE, 
normalizeTrace = TRUE, contrasts = NULL, 
model = TRUE, X = TRUE, Y = TRUE, K = TRUE, ...)
varComp.fit(Y, X = matrix(0, length(Y), 0L), K, control = varComp.control())
doFit.varComp(object)

Arguments

fixed

A two-sided formula specifying the fixed effects of the model. This is the same as in stats::lm or nlme::lme.

data

an optional data frame, list or environment containing the variables in the model. If not found in data, the variables are taken from environment(fixed), typically the environment from which varComp is called. This is the same as the data argument of stats::lm or nlme::lme.

random

An optional one-sided formula specifying additive random effects, of the form ~ z1 + z2 + z3. Interactions are allowed. The error variance component should not be included. A warning will be given when environment(fixed) and environment(random) do not match. See details.

varcov

An optional list of symmetric positive semidefinite matrices. If random is non-missing, these matrices represent the correlation matrix of each random effect. Thus the number of the random effects in random must be equal to the length of varcov. If varcov is named, the names will be matched to those used in random, following the same rule of arguments matching in function calls. If varcov is unnamed, it is assumed that the order is the same as in random. If random is missing, the weighted sum of these matrices represent the contribution of random effects to the marginal variance of the response variable, with unknown weights representing variance components.

weights

An optional nonnegative vector of the same length as the response variable specified in fixed. When it is given, it is inversely proportional to the error variances not captured by random and varcov. This is similar to the weights argument in stats::lm.

subset

An optional vector specifying a subset of observations to be used in the fitting process.

family

The same as the family argument of stats::glm. However, only gaussian('identity') is supported currently.

na.action

The same as in stats::lm.

offset

The same as in stats::glm. These offsets are assumed as known fixed effects.

control

An object from calling varComp.control.

doFit

A logical scalar, indicating whether model fitting should be performed.

normalizeTrace

A logical scalar, indicating whether the individual variance-covariance matrices should be normalized such that variance components are on the same scale.

contrasts

The same as in stats::lm.

model

A logical scalar, indicating whether the model frame will be included in the result.

For varComp, this is a logical scalar, indicating whether the fixed-effect design matrix should be included in the result. For varComp.fit this is the optional numeric fixed effect design matrix for the model. If X is missing or a matrix with zero columns, it is assumed that Y has zero mean.

For varComp, this is a logical scalar, indicating whether the response variable should be included in the result. For varComp.fit, this is a numeric vector of response variables.

For varComp, this is a logical scalar, indicating whether the list of variance-covariance matrices should be included in the result. These matrices are computed by pre-and-post multiplying the design matrices of random effects specified by random and the corresponding matrices specified by varcov. For varComp.fit, this is a list of variance-covariance matrices.

…

When control is given, this is ignored. Otherwise, these arguments are passed to varComp.control and the result will be used as the control argument.

object

An object of class varComp.

Value

The value of any of the three functions is a list with class varComp. Depending on whether doFit is TRUE or FALSE, the components in the list might differ. In either case, the result from varComp always contains the following elements:

na.action: the na.action used in the model frame.
offset: the same as input.
contrasts: the contrast used in the fixed-effect design matrix.
xzlevels: the levels of both fixed-effect and random-effect factors.
terms: both fixed-effect and random-effect terms.
call: the actual call.
nobs: the number of observations.
control: the same as input.
random.labels: the labels used to differentiate random effects. Note that the error variance is not included here. It is safe to check the number of variance components specified by the model by checking the length of random.labels.
doFit: The value of doFit is either TRUE or a call. If the doFit argument of varComp is set to FALSE, this component will be a call that will do the model fitting when being evaluated. If the model fitting has already been done, this component will always be TRUE.

In the case where doFit is TRUE, the following components will also be present in the result of varComp and doFit:

parms: a named numeric vector of parameter estimates. These are the estimated ratio of each variance component relative to the error variance.
gradients: a named numeric vector of gradient of the PREML function at the final parameter estimate. Because of non-negativity constraints, the gradients are not necessarily all zero at convergence.
hessian: a numeric matrix of Hessian matrix at parameter estimate. This is always computed through the observed information, no matter how the information argument of varComp.control is set.
sigma2: the estimated error variance.
varComps: a named numeric vector of variance components. These are the same as parms times sigma2.
n.iter: the number of iterations during fitting.
PREML: the final maximized PREML function value.
X.Q2: a matrix, when left-multiplied to the response variable, produces the residual contrasts.
residual.contrast: a numeric vector of residual contrasts with length{Y} - qr(X)$rank elements.
working.cor: the list of individual variance-covariance matrices, whose weighted sum is the marginal variance-covariance matrix of the residual contrast.

If varComp.fit is called directly, then doFit is always TRUE, but the result will not contain any of na.action, offset, contrasts, xzlevels, and terms, as they do not apply.

For the varComp interface, if any of model, X, Y or K is TRUE, the corresponding part will be included in the result. If weights is non-missing, it will be included in the result.

Details

The variance component model is of form $$\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{e}$$ where $\mathbf{e}$ is multivariate normally distributed with mean zero and variance-covariance matrix $\mathbf{V}$ being $$\mathbf{V} = \sum_{j=1}^R \sigma_j^2 \mathbf{K}_j + \sigma_e^2 \mathbf{W}$$ in which $\mathbf{K}_j$ are known positive semidefinite matrices and $\mathbf{W}$ is a known diagonal positive definite matrix. In the case of random-effect modeling, the $\mathbf{K}$ matrices are further given by $$\mathbf{K}_j = \mathbf{Z}_j \mathbf{G}_j \mathbf{Z}_j^\top$$ where $\mathbf{Z}_j$ is the design matrix for the $j$th random-effect factor and $\mathbf{G}_j$ the known variance-covariance matrix for the $j$th random-effect vector.

In the varComp formula interface, the $\mathbf{X}$ matrix and response variable are specified by the fixed argument. The optional random argument specifies the $\mathbf{Z}$ matrices. When random is missing, they are assumed to be identity matrices. The varcov argument specifies the $\mathbf{G}$ matrices. The weights argument specifies the $\mathbf{W}$ matrix.

Note that in random, kernel functions are allowed. For example, random = ~ ibs(SNP) is allowed to specify the similarities contributed by a matrix SNP through the identity-by-descent (IBS) kernel.

Fixed-by-random interactions are allowed in random. If F is a fixed factor with 2 levels, and F has already appeared in the fixed argument, then random = ~ ibs(SNP) + F:ibs(SNP) will specify two random effects. The first contributes to the marginal variance-covariance matrix through IBS(SNP); the other contributes to the marginal variance-covariance matrix through tcrossprod(X_F2) * IBS(SNP), where X_F2 is the 2nd column of the fixed-effect design matrix of factor F under the sum-to-zero contrast. Note that sum-to-zero contrast will be used even if another contrast is used for F in the fixed argument. This behavior corresponds to what we usually mean by fixed-by-random interaction. If one indeed needs to avoid sum-to-zero contrast in this case, it is only possible to do so by providing varcov directly without specifying random.

Additionally, the fixed-by-random interactions in random specified by the ``:'' symbol are actually treated as if they are specified by ``*''. In other words, random = ~ ibs(SNP) + F:ibs(SNP), random = ~ F:ibs(SNP), and random = ~ F*ibs(SNP) are actually treated as equivalent when F has already appeared on the right-side of fixed. Again this corresponds to the usual notion of fixed-by-random interaction. If such behavior is not wanted, one can directly give the intended varcov and leave random as missing.

Finally, intercept is not included in the random interface.

The model fitting process uses profiled restricted maximum likelihood (PREML), where the error variance is always profiled out of the REML likelihood. The non-negativity constraints of variance components are always imposed. See varComp.control for arguments controlling the modeling fitting.

References

Qu L, Guennel T, Marshall SL. (2013) Linear Score Tests for Variance Components in Linear Mixed Models and Applications to Genetic Association Studies. Biometrics, Volume 69, Issue 4, pages 883--892.

Examples

Run this code

# NOT RUN {
### Oxide/Semiconductor data example
library(nlme)
data(Oxide)
lmef = lme(Thickness~Source, Oxide, ~1|Lot/Wafer)
vcf = varComp(Thickness~Source, Oxide, ~Lot/Wafer)
VarCorr(lmef)
coef(vcf, 'varComp') ## same values as above
k0=tcrossprod(model.matrix(~0+Lot,Oxide))
k1=tcrossprod(Oxide$Source==1)*k0
k2=tcrossprod(Oxide$Source==2)*k0
k3=tcrossprod(model.matrix(~0+Lot:Wafer, Oxide))
## unequal variance across Source for Lot effects, in a preferred parameterization:
(vcf1 = varComp(Thickness~Source, Oxide, varcov=list(S1Lot=k1, S2Lot=k2, `Lot:Wafer`=k3))) 
## unequal variance across Source for Lot effects, in a different parameterization:
(vcf2 = varComp(Thickness~Source, Oxide, varcov=list(Lot=k0, S2Lot=k2, `Lot:Wafer`=k3)))
## unequal variance across Source for Lot effects, but in a poor parameterization that 
##   turns out to be the same as vcf after fitting.
(vcf3 = varComp(Thickness~Source, Oxide, varcov=list(Lot=k0, S1Lot=k1, `Lot:Wafer`=k3)))  
logLik(vcf)
logLik(vcf1)
logLik(vcf2)  ## the same as vcf1
logLik(vcf3)  ## the same as vcf
## fixef-effect only
vcf0 = varComp(Thickness~Source, Oxide)
summary(vcf0)
summary(lmf<-lm(Thickness~Source, Oxide))
vcf00 = varComp(Thickness~0, Oxide)
summary(vcf00)
summary(lmf0<-lm(Thickness~0, Oxide))



### Genetics example	
trt=gl(2, 15)
set.seed(2340)
dat=data.frame(trt=trt)
dat$SNP=matrix(sample(0:2, 120, replace=TRUE), 30)
dat$Y  = as.numeric(trt)+rnorm(30) + dat$SNP%*%rnorm(4)
(vcf0 = varComp(Y~trt, dat, ~ibs(SNP)))
(vcf00 = varComp(Y~trt, dat, varcov = list(`ibs(SNP)`=IBS(dat$SNP)))) ## same as above
(vcf1 = varComp(Y~trt, dat, ~ibs(SNP):trt)) ## two variance components

dat$trt[1]=NA
varComp(Y~trt, dat, ~ibs(SNP))  ## 29 observations compared to 30 in vcf0
dat$SNP[2,1]=NA
varComp(Y~trt, dat, ~ibs(SNP))  ## still, 29 observations, as ibs handles sporadic NA
dat$SNP[3, ]=NA
varComp(Y~trt, dat, ~ibs(SNP))  ## 28 observations, as no genotype for the 3rd obs



# }

Run the code above in your browser using DataLab