varComp: Fitting variance component models
In varComp: Variance Component Models

Description Usage Arguments Details Value Note Author(s) References See Also Examples

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.

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)

`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.
`X`	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.
`Y`	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.
`K`	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`.

The variance component model is of form

Y=X β + e

where e is multivariate normally distributed with mean zero and variance-covariance matrix V being

V= σ{j=1}^R σ_j^2 K_j + σ_e^2 W

in which K_j are known positive semidefinite matrices and W is a known diagonal positive definite matrix. In the case of random-effect modeling, the K matrices are further given by

K_j = Z_j G_j Z_j^T

where Z_j is the design matrix for the jth random-effect factor and G_j the known variance-covariance matrix for the jth random-effect vector.

In the varComp formula interface, the X matrix and response variable are specified by the fixed argument. The optional random argument specifies the Z matrices. When random is missing, they are assumed to be identity matrices. The varcov argument specifies the G matrices. The weights argument specifies the 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.

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.

Currently, available S3 methods for the varComp class are:

coef: retrieving fixed-effect estimates, variance component estimates, or estimates of the ratios of variance components to the error variance.
model.matrix: retrieving the fixed-effect design matrix, equivalent random-effect design matrix, or the list of individual variance-covariance matrices.
vcov: computing variance-covariance matrix of fixed-effect parameters, variance parameters, or the response variable.
fixef: retrieving fixed-effect parameter estimates with standard errors and default hypothesis tests against zero. See also satterth.varComp and KR.varComp for computing approximate denominator of degrees of freedom for the F-statistics.
anova: similar to fixef, but returning a data frame.
varComp.test: testing for nullity of random components (other than the error variance).
logLik: retrieving the maximized PREML value with degrees of freedom.
print: showing the object.
summary: summarize the object by providing standard errors and default tests against zero.

Long Qu

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.

varComp.test, coef.varComp, model.matrix.varComp, vcov.varComp, fixef.varComp, satterth.varComp, KR.varComp, logLik.varComp, print.varComp, summary.varComp, nlme::lme, stats::lm

### 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