Package 'VCA'

Description

This package implements ANOVA-type estimation of variance components (VC) for linear mixed models (LMM), and provides Restricted Maximum Likelihood (REML) estimation incorporating functionality of the lme4 package. For models fitted by REML the typical VCA-table is derived, also containing the variances of VC, which are approximated by the method outlined in Giesbrecht & Burns (1985). REML-estimation is available via functions remlVCA for variance component analysis (VCA) and remlMM for fitting general LMM.

ANOVA-methodology is a special method of moments approach for estimating (predicting) variance components implemented in functions anovaMM and anovaVCA. The former represents a general, unrestricted approach to fitting linear mixed models, whereas the latter is tailored for performing a VCA on random models. Experiments of this type frequently occur in performance evaluation analyses of diagnostic tests or analyzers (devices) quantifying various types of precision (see e.g. guideline EP05-A2/A3 of the Clinical and Laboratory Standards Institute - CLSI).

The general Satterthwaite approximation of denominator degrees of freedom for tests of fixed effects (test.fixef) and LS Means (test.lsmeans) is implemented as used in SAS PROC MIXED. Results differ for unbalanced designs because of the different approaches to estimating the covariance matrix of variance components. Here, two algorithms are implemented for models fitted via ANOVA, $1^{st}$ the "exact" method described in Searle et. al (1992), $2^{nd}$ an approximation described in Giesbrecht & Burns (1985). The latter is also used for models fitted by REML. See test.fixef and getGB for details on this topic.

Furthermore, the Satterthwaite approximation of degrees of freedom for individual VCs and total variance is implemented. These are employed in Chi-Squared tests of estimated variances against a claimed value (total, error), as well as in Chi-Squared based confidence intervals (CI) (see VCAinference). Whenever ANOVA-type estimated VCs become negative, the default is to set them equal to 0. ANOVA mean squares used within the Satterthwaite approximation will be adapted to this situation by re-computing ANOVA mean squares ($s_{MS}$) as $s_{MS} = C * \sigma^{2 }$, where $C$ is a coefficient matrix and a function of the design matrix and $\sigma^{2}$ is the column-vector of adapted variance components. Total variance corresponds to a conservative estimate of the total variability in these cases, i.e. it will be larger than e.g. the total variance of the same model fitted by REML, because the negative VC will not contribute to total variance. See the documentation anovaVCA and anovaMM for details, specifically argument NegVC.

Additionally to fitting linear mixed models and performing VCA-analyses, various plotting methods are implemented, e.g. a variability chart visualizing the variability in sub-classes emerging from an experimental design (varPlot). Random effects and residuals can be transformed and plotted using function plotRandVar. Standardization and studentization are generally available, Pearson-type transformation is only available for residuals. Plotting (studentized) random variates of a LMM should always be done to reveal potential problems of the fitted model, e.g. violation of model assumptions and/or whether there are outlying observations.

There are not any more two approaches to estimating ANOVA sums (SSQ) of squares as in previous package-versions. Now, only a fast FORTRAN-routine is used generating the column vector of SSQ, coefficient matrix C (previously computed using the Abbreviated Doolittle and Square Root Method), and the covariance matrix of VC. Overall, this leads to a dramatic reduction of computation time for models fitted using ANOVA.

Further reduction of the computation time can be achieved using Intel's Math Kernel Library (MKL). When the package is loaded it will be automatically checked whether this is the case or not.

In LS Means computation of fitted LMM it is possible to compute LS Means using specific values of covariables, which is equivalent to using option 'AT' in the 'lsmeans'-statement of SAS PROC MIXED. It is also possible to apply other than the default weighting scheme for (fixed) factor-variables. See the details section in lsmeans and the description of argument at.

Note: The 'UnitTests' directory within the package-directory contains a pre-defined test-suite which can be run by sourcing 'RunAllTests.R' for user side testing (installation verification). It requires the 'RUnit' package and checks the numerical equivalence to reference results (SAS PROC MIXED method=type1/reml, SAS PROC VARCOMP) for balanced and unbalanced data and different experimental designs.

Details

Author(s)

Andre Schuetzenmeister [email protected], Florian Dufey [email protected]

References

Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York

Goodnight, J.H. (1979), A Tutorial on the SWEEP Operator, The American Statistician, 33:3, 149-158

Giesbrecht, F.G. and Burns, J.C. (1985), Two-Stage Analysis Based on a Mixed Model: Large-Sample Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486

Satterthwaite, F.E. (1946), An Approximate Distribution of Estimates of Variance Components., Biometrics Bulletin 2, 110-114

Gaylor,D.W., Lucas,H.L., Anderson,R.L. (1970), Calculation of Expected Mean Squares by the Abbreviated Doolittle and Square Root Methods., Biometrics 26 (4): 641-655

SAS Help and Documentation PROC MIXED, SAS Institute Inc., Cary, NC, USA

Description

Estimate/Predict random effects employing ANOVA-type estimation and obtain generalized least squares estimates of fixed effects for any linear mixed model including random models and linear models.

Usage

anovaMM(
  form,
  Data,
  by = NULL,
  VarVC.method = c("scm", "gb"),
  NegVC = FALSE,
  quiet = FALSE,
  order.data = TRUE
)

Arguments

Details

A Linear Mixed Model, noted in standard matrix notation, can be written as $y = Xb + Zg + e$, where $y$ is the column vector of observations, $X$ and $Z$ are design matrices assigning fixed ($b$), respectively, random ($g$) effects to observations, and $e$ is the column vector of residual errors. Whenever there is an intercept in the model, i.e. the substring "-1" is not part of the model formula, the same restriction as in SAS PROC MIXED is introduced setting the last fixed effect equal to zero. Note, that the results of an linear contrasts are not affected by using an intercept or not, except that constrained fixed effects cannot be part of such contrasts (one could use the intercept estimated instead).

Here, no further restrictions on the type of model are made. One can fit mixed models as well as random models, which constitute a sub-set of mixed models (intercept being the only fixed effect). Variables must be either of type "numeric" or "factor". "character" variables are automatically converted to factors and the response variable has to be numeric, of course. In case that 'class(Data[,i])' is neither one of these three options, an error is issued. Even simple linear models can be fitted, i.e. models without a random part (without $Zg$) besides the residual errors. In this case, an Analysis of Variance (ANOVA) table is computed in the same way as done by function 'anova.lm'.

One drawback of using ANOVA-type estimation of random effects is, that random effects are independent, i.e they have zero covariance by definition $cov(g_i,g_j) = 0$. Another one is that estimated variance components may become negative under certain conditions. The latter situation is addressed by setting negative variance estimates equal to zero and adapting ANOVA mean squares (MS) as $MS = C * VC$, where $C$ is a coefficient matrix and a function of the design matrix $[X Z]$ and $VC$ is the column-vector of adapted variance components. The Satterthwaite approximation of total degrees of freedom (DF for total variance) will use adapted $MS$-values.

Note, that setting negative VCs equal to zero results in a conservative estimate of the total variance, i.e. it will be larger than the estimate including the negative VC(s). Use parameter 'NegVC=TRUE' to explicitly allow negative variance estimates.

For further details on ANOVA Type-I estimation methods see anovaVCA.

Value

(VCA) object

Author(s)

Andre Schuetzenmeister [email protected]

References

Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York

Goodnight, J.H. (1979), A Tutorial on the SWEEP Operator, The American Statistician, 33:3, 149-158

Giesbrecht, F.G. and Burns, J.C. (1985), Two-Stage Analysis Based on a Mixed Model: Large-Sample Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486

H.P.Piepho, A.Buechse and K.Emrich (2003), A Hitchhiker's Guide to Mixed Models for Randomized Experiments, J.Agronomy & Crop Science 189, p. 310-322

Gaylor,D.W., Lucas,H.L., Anderson,R.L. (1970), Calculation of Expected Mean Squares by the Abbreviated Doolittle and Square Root Methods., Biometrics 26 (4): 641-655

SAS Help and Documentation PROC MIXED, SAS Institute Inc., Cary, NC, USA

See Also

anovaVCA

anovaVCA, VCAinference, remlVCA, remlMM ranef, fixef, vcov, vcovVC, test.fixef, test.lsmeans, plotRandVar

Examples

## Not run: 

data(dataEP05A2_2)

# assuming 'day' as fixed, 'run' as random
anovaMM(y~day/(run), dataEP05A2_2)

# assuming both as random leads to same results as
# calling anovaVCA
anovaMM(y~(day)/(run), dataEP05A2_2)
anovaVCA(y~day/run, dataEP05A2_2)

# use different approaches to estimating the covariance of 
# variance components (covariance parameters)
dat.ub <- dataEP05A2_2[-c(11,12,23,32,40,41,42),]			# get unbalanced data
m1.ub <- anovaMM(y~day/(run), dat.ub, VarVC.method="scm")
m2.ub <- anovaMM(y~day/(run), dat.ub, VarVC.method="gb")		
V1.ub <- round(vcovVC(m1.ub), 12)
V2.ub <- round(vcovVC(m2.ub), 12)
all(V1.ub == V2.ub)

# fit a larger random model
data(VCAdata1)
fitMM1 <- anovaMM(y~((lot)+(device))/(day)/(run), VCAdata1[VCAdata1$sample==1,])
fitMM1
# now use function tailored for random models
fitRM1 <- anovaVCA(y~(lot+device)/day/run, VCAdata1[VCAdata1$sample==1,])
fitRM1

# there are only 3 lots, take 'lot' as fixed 
fitMM2 <- anovaMM(y~(lot+(device))/(day)/(run), VCAdata1[VCAdata1$sample==2,])

# the following model definition is equivalent to the one above,
# since a single random term in an interaction makes the interaction
# random (see the 3rd reference for details on this topic)
fitMM3 <- anovaMM(y~(lot+(device))/day/run, VCAdata1[VCAdata1$sample==2,])

# fit same model for each sample using by-processing
lst <- anovaMM(y~(lot+(device))/day/run, VCAdata1, by="sample")
lst

# fit mixed model originally from 'nlme' package

library(nlme)
data(Orthodont)
fit.lme <- lme(distance~Sex*I(age-11), random=~I(age-11)|Subject, Orthodont) 

# re-organize data for using 'anovaMM'
Ortho <- Orthodont
Ortho$age2 <- Ortho$age - 11
Ortho$Subject <- factor(as.character(Ortho$Subject))
fit.anovaMM1 <- anovaMM(distance~Sex*age2+(Subject)*age2, Ortho)

# use simplified formula avoiding unnecessary terms
fit.anovaMM2 <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho)

# and exclude intercept
fit.anovaMM3 <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho)

# compare results
fit.lme
fit.anovaMM1
fit.anovaMM2
fit.anovaMM3

# are there a sex-specific differences?
cmat <- getL(fit.anovaMM3, c("SexMale-SexFemale", "SexMale:age2-SexFemale:age2")) 
cmat

test.fixef(fit.anovaMM3, L=cmat)

# former versions of the package used R-function 'lm' and 'anova',
# which is significantly slower for sufficiently large/complex models
data(realData)
datP1 <- realData[realData$PID==1,]
system.time(anova.lm.Tab <- anova(lm(y~lot/calibration/day/run, datP1)))
# Using the sweeping approach for estimating ANOVA Type-1 sums of squares
# this is now the default setting. 
system.time(anovaMM.Tab1  <- anovaMM(y~lot/calibration/day/run, datP1))

# compare results, note that the latter corresponds to a linear model,
# i.e. without random effects. Various matrices have already been computed,
# e.g. "R", "V" (which are identical in this case).
anova.lm.Tab
anovaMM.Tab1

## End(Not run)

Description

This function equates observed ANOVA Type-I sums of squares ($SS$) to their expected values and solves the resulting system of linear equations for variance components.

Usage

anovaVCA(
  form,
  Data,
  by = NULL,
  NegVC = FALSE,
  VarVC.method = c("scm", "gb"),
  MME = FALSE,
  quiet = FALSE,
  order.data = TRUE
)

Arguments

Details

For diagnostics, a key parameter is "precision", i.e. the accuracy of a quantification method influenced by varying sources of random error. This type of experiments is requested by regulatory authorities to proof the quality of diagnostic tests, e.g. quantifying intermediate precision according to CLSI guideline EP5-A2/A3. No, fixed effects are allowed besides the intercept. Whenever fixed effects are part of the model to be analyzed, use function anovaMM instead.

Function anovaVCA is tailored for performing Variance Component Analyses (VCA) for random models, assuming all VCs as factor variables, i.e. their levels correspond to distinct columns in the design matrix (dummy variables). Any predictor variables are automatically converted to factor variables, since continuous variables may not be used on the right side of the formula 'form'.

ANOVA $SS$ are computed employing the SWEEP-operator (Goodnight 1979, default). according to Searle et al. (1992) which corresponds to VarVC.method="scm".

Function anovaVCA represents a special form of the "method of moments" approach applicable to arbitrary random models either balanced or unbalanced. The system of linear equations, which is built from the ANOVA Type-I sums of squares, is closely related to the method used by SAS PROC VARCOMP, where ANOVA mean squares ($MS$) are used. The former can be written as $ss = C * s$ and the latter as $ms = D * s$, where $C$ and $D$ denote the respective coefficient matrices, $s$ the column-vector of variance components (VC) to be estimated/predicted, and $ss$ and $ms$ the column vector of ANOVA sum of squares, respectively, mean squares. Mutliplying element $d_ij$ of matrix $D$ by element $c_in$ of matrix $C$ ($i,j = 1,...,n$), results in matrix $C$. Thus, $C$ can easily be converted to $D$ by the inverse operation. Matrix $D$ is used to estimate total degrees of freedom (DF) according to Satterthwaite (1946).

The method for computing ANOVA Type-I $SS$ is much faster than fitting the linear model via lm and calling function anova on the 'lm' object for complex models, where complex refers to the number of columns of the design matrix and the degree of unbalancedness. $DF$ are directly derived from the SWEEP-operator as the number of linearly independent columns of the partial design matrix corresponding to a specific $VC$.

Value

(object) of class 'VCA'

Author(s)

Andre Schuetzenmeister [email protected]

References

Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York

Goodnight, J.H. (1979), A Tutorial on the SWEEP Operator, The American Statistician, 33:3, 149-158

Giesbrecht, F.G. and Burns, J.C. (1985), Two-Stage Analysis Based on a Mixed Model: Large-Sample Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486

Satterthwaite, F.E. (1946), An Approximate Distribution of Estimates of Variance Components., Biometrics Bulletin 2, 110-114

Gaylor,D.W., Lucas,H.L., Anderson,R.L. (1970), Calculation of Expected Mean Squares by the Abbreviated Doolittle and Square Root Methods., Biometrics 26 (4): 641-655

SAS Help and Documentation PROC MIXED, SAS Institute Inc., Cary, NC, USA

SAS Help and Documentation PROC VARCOMP, SAS Institute Inc., Cary, NC, USA

See Also

anovaMM, remlVCA, remlMM, print.VCA, VCAinference, ranef, plotRandVar, stepwiseVCA

Examples

## Not run: 

# load data (CLSI EP05-A2 Within-Lab Precision Experiment) 
data(dataEP05A2_2)

# perform ANOVA-estimation of variance components
res <- anovaVCA(y~day/run, dataEP05A2_2)
res

# design with two main effects (ignoring the hierarchical structure of the design)
anovaVCA(y~day+run, dataEP05A2_2)

# compute confidence intervals, perform F- and Chi-Squared tests
INF <- VCAinference(res, total.claim=3.5, error.claim=2)
INF

### load data from package
data(VCAdata1)

data_sample1 <- VCAdata1[VCAdata1$sample==1,]

### plot data for visual inspection
varPlot(y~lot/day/run, data_sample1)

### estimate VCs for 4-level hierarchical design (error counted) for sample_1 data
anovaVCA(y~lot/day/run, data_sample1)

### using different model (ignoring the hierarchical structure of the design)
anovaVCA(y~lot+day+lot:day:run, data_sample1)

### same model with unbalanced data
anovaVCA(y~lot+day+lot:day:run, data_sample1[-c(1,11,15),])

### use the numerical example from the CLSI EP05-A2 guideline (p.25)
data(Glucose,package="VCA")
res.ex <- anovaVCA(result~day/run, Glucose)

### also perform Chi-Squared tests
### Note: in guideline claimed SD-values are used, here, claimed variances are used
VCAinference(res.ex, total.claim=3.4^2, error.claim=2.5^2)

### now use the six sample reproducibility data from CLSI EP5-A3
### and fit per sample reproducibility model
data(CA19_9)
fit.all <- anovaVCA(result~site/day, CA19_9, by="sample")

reproMat <- data.frame(
Sample=c("P1", "P2", "Q3", "Q4", "P5", "Q6"),
Mean= c(fit.all[[1]]$Mean, fit.all[[2]]$Mean, fit.all[[3]]$Mean, 
fit.all[[4]]$Mean, fit.all[[5]]$Mean, fit.all[[6]]$Mean),
Rep_SD=c(fit.all[[1]]$aov.tab["error","SD"], fit.all[[2]]$aov.tab["error","SD"],
fit.all[[3]]$aov.tab["error","SD"], fit.all[[4]]$aov.tab["error","SD"],
fit.all[[5]]$aov.tab["error","SD"], fit.all[[6]]$aov.tab["error","SD"]),
Rep_CV=c(fit.all[[1]]$aov.tab["error","CV[%]"],fit.all[[2]]$aov.tab["error","CV[%]"],
fit.all[[3]]$aov.tab["error","CV[%]"],fit.all[[4]]$aov.tab["error","CV[%]"],
fit.all[[5]]$aov.tab["error","CV[%]"],fit.all[[6]]$aov.tab["error","CV[%]"]),
WLP_SD=c(sqrt(sum(fit.all[[1]]$aov.tab[3:4,"VC"])),sqrt(sum(fit.all[[2]]$aov.tab[3:4, "VC"])),
sqrt(sum(fit.all[[3]]$aov.tab[3:4,"VC"])),sqrt(sum(fit.all[[4]]$aov.tab[3:4, "VC"])),
sqrt(sum(fit.all[[5]]$aov.tab[3:4,"VC"])),sqrt(sum(fit.all[[6]]$aov.tab[3:4, "VC"]))),
WLP_CV=c(sqrt(sum(fit.all[[1]]$aov.tab[3:4,"VC"]))/fit.all[[1]]$Mean*100,
sqrt(sum(fit.all[[2]]$aov.tab[3:4,"VC"]))/fit.all[[2]]$Mean*100,
sqrt(sum(fit.all[[3]]$aov.tab[3:4,"VC"]))/fit.all[[3]]$Mean*100,
sqrt(sum(fit.all[[4]]$aov.tab[3:4,"VC"]))/fit.all[[4]]$Mean*100,
sqrt(sum(fit.all[[5]]$aov.tab[3:4,"VC"]))/fit.all[[5]]$Mean*100,
sqrt(sum(fit.all[[6]]$aov.tab[3:4,"VC"]))/fit.all[[6]]$Mean*100),
Repro_SD=c(fit.all[[1]]$aov.tab["total","SD"],fit.all[[2]]$aov.tab["total","SD"],
fit.all[[3]]$aov.tab["total","SD"],fit.all[[4]]$aov.tab["total","SD"],
fit.all[[5]]$aov.tab["total","SD"],fit.all[[6]]$aov.tab["total","SD"]),
Repro_CV=c(fit.all[[1]]$aov.tab["total","CV[%]"],fit.all[[2]]$aov.tab["total","CV[%]"],
fit.all[[3]]$aov.tab["total","CV[%]"],fit.all[[4]]$aov.tab["total","CV[%]"],
fit.all[[5]]$aov.tab["total","CV[%]"],fit.all[[6]]$aov.tab["total","CV[%]"]))

for(i in 3:8) reproMat[,i] <- round(reproMat[,i],digits=ifelse(i%%2==0,1,3))
reproMat

# now plot the precision profile over all samples
plot(reproMat[,"Mean"], reproMat[,"Rep_CV"], type="l", main="Precision Profile CA19-9",
xlab="Mean CA19-9 Value", ylab="CV[%]")
grid()
points(reproMat[,"Mean"], reproMat[,"Rep_CV"], pch=16)

# load another example dataset and extract the "sample==1" subset
data(VCAdata1)
sample1 <- VCAdata1[which(VCAdata1$sample==1),]

# generate an additional factor variable and random errors according to its levels
sample1$device <- gl(3,28,252)                                      
set.seed(505)
sample1$y <- sample1$y + rep(rep(rnorm(3,,.25), c(28,28,28)),3)     

# fit a crossed-nested design with main factors 'lot' and 'device' 
# and nested factors 'day' and 'run' nested below 
res1 <- anovaVCA(y~(lot+device)/day/run, sample1) 
res1

# fit same model for each sample using by-processing
lst <- anovaVCA(y~(lot+device)/day/run, VCAdata1, by="sample")
lst

# now fitting a nonsense model on the complete dataset "VCAdata1" 
# where the SWEEP-operator is the new default since package version 1.2
# takes ~5s
system.time(res.sw <- anovaVCA(y~(sample+lot+device)/day/run, VCAdata1))
# applying functions 'anova' and 'lm' in the same manner takes ~ 265s
system.time(res.lm <- anova(lm(y~(sample+lot+device)/day/run, VCAdata1)))
res.sw
res.lm

## End(Not run)

Description

Standard 'as.matrix' Method for 'VCA' S3-Objects

Usage

## S3 method for class 'VCA'
as.matrix(x, ...)

Arguments

Value

(matrix) equal to x$aov.tab with additional attributes "Mean" and "Nobs"

Author(s)

Andre Schuetzenmeister [email protected]

See Also

as.matrix.VCAinference

Examples

## Not run: 
data(dataEP05A2_1)
fit <- anovaVCA(y~day/run, dataEP05A2_1)
as.matrix(fit)

## End(Not run)

Description

This function makes use of the hidden feature of function print.VCAinference which invisibly returns character matrices of estimated variance components expressed as "VC" (variance component), "SD" (standard deviation) or "CV" (coefficient of variation). If argument "what" is not specified, a named list will be returned with all three matrices.

Usage

## S3 method for class 'VCAinference'
as.matrix(x, what = c("VC", "SD", "CV"), digits = 6, ...)

Arguments

Value

(matrix) with point estimates, one- and two-sided confidence intervals and variances of the estimated variance components

Author(s)

Andre Schuetzenmeister [email protected]

See Also

print.VCAinference, as.matrix.VCA

Examples

## Not run: 
data(dataEP05A2_1)
fit <- anovaVCA(y~day/run, dataEP05A2_1)
inf <- VCAinference(fit, VarVC=TRUE)
as.matrix(inf, what="VC", digits=6)
as.matrix(inf, what="SD", digits=6)
as.matrix(inf, what="CV", digits=2)

# request list of matrices
as.matrix(inf)

## End(Not run)

Description

Function buildList creates a nested-list reflecting the hierarchical structure of a fully-nested model, respectively, the imputed hierarchical structure of the data (see details).

Usage

buildList(
  Data,
  Nesting,
  Current,
  resp,
  keep.order = TRUE,
  useVarNam = TRUE,
  sep = "",
  na.rm = TRUE,
  Points = list(pch = 16, cex = 0.5, col = "black")
)

Arguments

Details

This function is not intended to be used directly and serves as helper function for varPlot. Each factor-level, on each level of nesting is accompanied with a set of descriptive statistics, such as mean, median, var, sd, ... which can be evaluated later on. These information are used in function varPlot, which implements a variability chart. Note, that this function is also used if data does not correspond to a fully-nested design, i.e. the hierarchical structure is inferred from the model formula. The order of main factors (not nested within other factors) appearing in the model formula determines the nesting structure imputed in order to plot the data as variability chart.

Value

(list) which was recursively built, representing the data of the fully-nested as hierarchy

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 
# load data (CLSI EP05-A2 Within-Lab Precision Experiment)
data(dataEP05A2_3)

# build a list representing the hierarichal structure of a fully-nested model
# there needs to be a distinct hierarchy for being able to plot the data
# as variability chart (this function is not exported)
lst <- VCA:::buildList(Data=dataEP05A2_3, Nesting=c("day", "run"), Current="day", resp="y")

## End(Not run)

Description

This data set consists of the example data of a complete three-site, multisample reproducibility study as presented in the CLSI EP5-A3 guideline. It shows quantitative results of an automated immunometric assay measuring parameter CA19-9. This dataset is described in Appendix B of this guideline consisting of 6 samples, each measured on one of three sites, at five days with five replicates per day.

Usage

data(CA19_9)

Format

data.frame with 450 rows and 4 variables.

References

CLSI EP05-A3 - Evaluation of Precision of Quantitative Measurement Procedures; Approved Guideline - Third Edition. CLSI

Description

Majority of the code is borrowed from the Microsoft R Open Rprofile.site file. In case MKL can be detected this information will be stored in a separate envrionment, which other function know about. If so, an optimized version of function getGB will be used which used ordinary matrix-objects instead of matrices defined by the Matrix-package. This seems to accelerate computation time for large datasets by up to factor 30.

Usage

check4MKL()

Details

This function is for internal use only and therefore not exported.

Value

variable 'MKL' in envir "msgEnv" will be set to TRUE/FALSE

Author(s)

Authors of the Rprofile.site file in Microsoft R Open, Andre Schuetzenmeister [email protected]

Description

This function is intended to check a variance component analysis either before or after performing it. This is particularily important for less experienced users who my not exactly know where error messages come from. External software using functions anovaVCA or remlVCA also via function fitVCA may also benefit from more user-friendly error messages.

Usage

checkData(form, Data)

Arguments

Value

(list) of length equal to the number of terms in 'form' each element being a list of messages with identified, potential problems.

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 
data(dataEP05A2_1)
dat0 <- dataEP05A2_1[1:16,]
# everything should be ok
checkData(y~day/run, dat0)
# data identical response for all obs 
dat1 <- dat0
dat1$y <- dat1[1,"y"]
remlVCA(y~day/run, dat1)
checkData(y~day/run, dat1)
# now factor-levels have identical values
dat2 <- dat0
dat2$y <- dat2$y[rep(seq(1,7,2), rep(2,4))] 
checkData(y~day/run, dat2)
remlVCA(y~day/run, dat2, quiet=TRUE)
# indistinguishable terms are also problematic
dat3 <- data.frame(	y=rnorm(8,10),
						day=paste("day",rep(c(1,2),c(4,4))), 
						run=rep(c(2,1), c(4,4)))
checkData(y~day/run, dat3)
anovaVCA(y~day/run, dat3)
# no replicates, thus, no error variability
dat4 <- dat0[seq(1,15,2),]
dat4$day <- factor(dat4$day)
dat4$run <- factor(dat4$run)
checkData(y~day/run, dat4)
remlVCA(y~day/run, dat4)

## End(Not run)

Description

Is is checked whether 'term2' is different from 'term1' in adding information to the model. If both are main factors, i.e. no interactions terms, it is checked whether levels of 'term2' differ from those of 'term1'. Otherwise, 'term2' is an interaction with a part being different from 'term1'.

Usage

checkVars(Data, term1, term2)

Arguments

Value

(list) with components \"Diff\"=part of 'term2' distinguishing it from 'term1', \"AddInfo\"=message informing about potential problems with 'term2'

Author(s)

Andre Schuetzenmeister [email protected]

Description

This dataset was added because it generates an error in function 'chol2inv' when trying to invert the variance-covariance matrix 'V' of the mixed model 'value~ID+(Site)'. This dataset and the associated mixed model are part of the unit-test collection of the package.

Usage

data(chol2invData)

Format

A data frame with 158 observations on the following 3 variables.

• value
  
  The response variable.
  
  

• ID
  
  Variable with 6 levels corresponding to samples.
  
  

• Site
  
  Variable with 3 levels corresponding to sites/devices.
  

Description

For conveniently using objects of class 'VCA' with other packages expecting this function, e.g. the 'multcomp' package for general linear hypotheses for parametric models (currently not fully implemented).

Usage

## S3 method for class 'VCA'
coef(object, quiet = FALSE, ...)

Arguments

Examples

## Not run: 
data(dataEP05A2_1)
fit1 <- anovaMM(y~day/(run), dataEP05A2_1)
coef(fit1)
fit2 <- anovaVCA(y~day/run, dataEP05A2_1)
coef(fit2)

## End(Not run)

Description

This data set consists of simulated measurements for an experiment conducted to evaluate the precision performance of measurement methods. On 20 days two separate runs with two replicates of the same sample are measured. Thus, factor 'day' is the top-level random factor (variance component), factor 'run' is nested within 'day'.

Usage

data(dataEP05A2_1)

Format

data.frame with 80 rows and 3 variables.

References

CLSI EP05-A2 - Evaluation of Precision Performance of Quantitative Measurement Methods. CLSI

Description

This data set consists of simulated measurements for an experiment conducted to evaluate the precision performance of measurement methods. On 20 days two separate runs with two replicates of the same sample are measured. Thus, factor 'day' is the top-level random factor (variance component), factor 'run' is nested within 'day'.

Usage

data(dataEP05A2_2)

Format

data.frame with 80 rows and 3 variables.

References

CLSI EP05-A2 - Evaluation of Precision Performance of Quantitative Measurement Methods. CLSI

Description

This data set consists of simulated measurements for an experiment conducted to evaluate the precision performance of measurement methods. On 20 days two separate runs with two replicates of the same sample are measured. Thus, factor 'day' is the top-level random factor (variance component), factor 'run' is nested within 'day'.

Usage

data(dataEP05A2_3)

Format

data.frame with 80 rows and 3 variables.

References

CLSI EP05-A2 - Evaluation of Precision Performance of Quantitative Measurement Methods. CLSI

Description

This data set consists of simulated measurements for an experiment to be conducted for evaluation of the precision performance of measurement methods. On 3 sites, on 5 days 5 replicates of the same sample are measured. Thus, factor 'site' is the top-level random factor (variance component), factor 'day' is nested within 'site'.

Usage

data(dataEP05A3_MS_1)

Format

data.frame with 75 rows and 3 variables.

References

Draft of CLSI EP05-A3 - Evaluation of Precision Performance of Quantitative Measurement Methods. CLSI

Description

This data set consists of simulated measurements for an experiment to be conducted for evaluation of the precision performance of measurement methods. On 3 sites, on 5 days 5 replicates of the same sample are measured. Thus, factor 'site' is the top-level random factor (variance component), factor 'day' is nested within 'site'.

Usage

data(dataEP05A3_MS_2)

Format

data.frame with 75 rows and 3 variables.

References

Draft of CLSI EP05-A3 - Evaluation of Precision Performance of Quantitative Measurement Methods. CLSI

Description

This data set consists of simulated measurements for an experiment to be conducted for evaluation of the precision performance of measurement methods. On 3 sites, on 5 days 5 replicates of the same sample are measured. Thus, factor 'site' is the top-level random factor (variance component), factor 'day' is nested within 'site'.

Usage

data(dataEP05A3_MS_3)

Format

data.frame with 75 rows and 3 variables.

References

Draft of CLSI EP05-A3 - Evaluation of Precision Performance of Quantitative Measurement Methods. CLSI

Description

This data set consists of 21 measurements of the same sample, suitable to quantify the mesurement error on the same device.

Usage

data(dataRS0003_1)

Format

data.frame with 21 rows and 1 variable.

Description

This data set consists of 21 measurements of the same sample, suitable to quantify the mesurement error on the same device.

Usage

data(dataRS0003_2)

Format

data.frame with 21 rows and 1 variable.

Description

This data set consists of 21 measurements of the same sample, suitable to quantify the mesurement error on the same device.

Usage

data(dataRS0003_3)

Format

data.frame with 21 rows and 1 variable.

Description

This data set consists of 15 measurements of the same sample, measured on 5 days with 3 measurements on each day. This small experiment is suitable to quantify between-day variability on the same device.

Usage

data(dataRS0005_1)

Format

data.frame with 15 rows and 2 variables.

Description

This data set consists of 15 measurements of the same sample, measured on 5 days with 3 measurements on each day. This small experiment is suitable to quantify between-day variability on the same device.

Usage

data(dataRS0005_2)

Format

data.frame with 15 rows and 2 variables.

Description

This data set consists of 15 measurements of the same sample, measured on 5 days with 3 measurements on each day. This small experiment is suitable to quantify between-day variability on the same device.

Usage

data(dataRS0005_3)

Format

data.frame with 15 rows and 2 variables.

Description

This is a helper function for function test.fixef approximating degrees of freedom for linear contrasts of fixed effect parameter estimates.

Usage

DfSattHelper(obj, x)

Arguments

Value

(matrix) corresponding to the variance-covariance matrix of fixed effects

Description

Function takes one or multiple objects and converts them to a single error-message. Objects can be output of functions try or checkData.

Usage

errorMessage(...)

Arguments

Value

(characer) string combining information from input-objects

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 
data(dataEP05A2_1)
dat2 <- dataEP05A2_1[1:16,]
dat2$y <- dat2$y[rep(seq(1,7,2), rep(2,4))] 
errorMessage(try(1/"a"), checkData(y~day/run, dat2))

## End(Not run)

Description

Function serves as interface to functions anovaMM and remlMM for fitting a linear mixed model (LMM) either by ANOVA or REML. All arguments applicable to either one of these functions can be specified (see anovaMM or remlMM for details).

Usage

fitLMM(
  form,
  Data,
  method = c("anova", "reml"),
  scale = TRUE,
  VarVC = TRUE,
  ...
)

Arguments

Details

Besides offering a convenient interface to both functions for fitting a LMM, this function also provides all elements required for standard task of fitted models, e.g. prediction, testing general linear hypotheses via R-package multcomp, etc. (see examples).

Author(s)

Andre Schuetzenmeister [email protected]

See Also

fitVCA, anovaMM, remlMM

Examples

## Not run: 
data(dataEP05A2_2)

# assuming 'day' as fixed, 'run' as random
# Note: default method is "anova"
fitLMM(y~day/(run), dataEP05A2_2)

# explicitly request "reml"
fitLMM(y~day/(run), dataEP05A2_2, method="reml")

# assuming both as random leads to same results as
# calling anovaVCA (ANOVA is the default)
fitLMM(y~(day)/(run), dataEP05A2_2)
anovaVCA(y~day/run, dataEP05A2_2)

# now using REML-estimation
fitLMM(y~(day)/(run), dataEP05A2_2, "reml")
remlVCA(y~day/run, dataEP05A2_2)

# use different approaches to estimating the covariance of 
# variance components (covariance parameters)
# create unbalanced data
dat.ub <- dataEP05A2_2[-c(11,12,23,32,40,41,42),]
m1.ub <- fitLMM(y~day/(run), dat.ub, VarVC.method="scm")
# VarVC.method="gb" is an approximation not relying on quadratic forms
m2.ub <- fitLMM(y~day/(run), dat.ub, VarVC.method="gb")		
# REML-estimated variance components usually differ from ANOVA-estimates
# and so do the variance-covariance matrices
m3.ub <- fitLMM(y~day/(run), dat.ub, "reml", VarVC=TRUE)		
V1.ub <- round(vcovVC(m1.ub), 12)
V2.ub <- round(vcovVC(m2.ub), 12)
V3.ub <- round(vcovVC(m3.ub), 12)

# fit a larger random model
data(VCAdata1)
fitMM1 <- fitLMM(y~((lot)+(device))/(day)/(run), VCAdata1[VCAdata1$sample==1,])
fitMM1
# now use function tailored for random models
fitRM1 <- anovaVCA(y~(lot+device)/day/run, VCAdata1[VCAdata1$sample==1,])
fitRM1

# there are only 3 lots, take 'lot' as fixed 
fitMM2 <- fitLMM(y~(lot+(device))/(day)/(run), VCAdata1[VCAdata1$sample==2,])
# use REML on this (balanced) data
fitMM2.2 <- fitLMM(y~(lot+(device))/(day)/(run), VCAdata1[VCAdata1$sample==2,], "reml")

# the following model definition is equivalent to the one above,
# since a single random term in an interaction makes the interaction
# random (see the 3rd reference for details on this topic)
fitMM3 <- fitLMM(y~(lot+(device))/day/run, VCAdata1[VCAdata1$sample==2,])

# fit same model for each sample using by-processing
lst <- fitLMM(y~(lot+(device))/day/run, VCAdata1, by="sample")
lst

# fit mixed model originally from 'nlme' package

library(nlme)
data(Orthodont)
fit.lme <- lme(distance~Sex*I(age-11), random=~I(age-11)|Subject, Orthodont) 

# re-organize data
Ortho <- Orthodont
Ortho$age2 <- Ortho$age - 11
Ortho$Subject <- factor(as.character(Ortho$Subject))
fit.anovaMM1 <- fitLMM(distance~Sex*age2+(Subject)*age2, Ortho)

# use simplified formula avoiding unnecessary terms
fit.anovaMM2 <- fitLMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho)

# and exclude intercept
fit.anovaMM3 <- fitLMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho)

# compare results
fit.lme
fit.anovaMM1
fit.anovaMM2
fit.anovaMM3

# are there a sex-specific differences?
cmat <- getL(fit.anovaMM3, c("SexMale-SexFemale", "SexMale:age2-SexFemale:age2")) 
cmat

test.fixef(fit.anovaMM3, L=cmat)

# fit LMM with fixed lot and device effects and test for lot-differences
data(VCAdata1)
fitS5 <- fitLMM(y~(lot+device)/(day)/(run), subset(VCAdata1, sample==5), "reml")
fitS5

# apply Tukey-HSD test to screen for lot differences
library(multcomp)
res.tuk <- glht(fitS5, linfct=mcp(lot="Tukey"))
summary(res.tuk)

# compact letter display
res.tuk.cld <- cld(res.tuk, col=paste0("gray", c(90,60,75)))
plot(res.tuk.cld)

## End(Not run)

Description

Function serves as interface to functions anovaVCA and remlVCA for fitting a variance component models (random models) either by ANOVA or REML. All arguments applicable to either one of these functions can be specified (see anovaVCA or remlVCA for details).

Usage

fitVCA(
  form,
  Data,
  method = c("anova", "reml"),
  scale = TRUE,
  VarVC = TRUE,
  ...
)

Arguments

Author(s)

Andre Schuetzenmeister [email protected]

See Also

fitLMM, anovaVCA, remlVCA

Examples

## Not run: 
#load data (CLSI EP05-A2 Within-Lab Precision Experiment) 
data(dataEP05A2_2)

# perform ANOVA-estimation of variance components
res.anova <- fitVCA(y~day/run, dataEP05A2_2, "anova")
# perform REML-estimation of variance components
res.reml <- fitVCA(y~day/run, dataEP05A2_2, "reml")

# compare scaling vs. not scaling the response
fit0 <- fitVCA(y~day/run, dataEP05A2_2, "anova", scale=TRUE)
fit1 <- fitVCA(y~day/run, dataEP05A2_2, "anova", scale=FALSE)

## End(Not run)

Description

Generic Method for Extracting Fixed Effects from a Fitted Model

Usage

fixef(object, ...)

Arguments

See Also

fixef.VCA

Description

Conveniently extracting the 'FixedEffects' element of an 'VCA' object.

Usage

## S3 method for class 'VCA'
fixef(
  object,
  type = c("simple", "complex"),
  ddfm = c("contain", "residual", "satterthwaite"),
  tol = 1e-12,
  quiet = FALSE,
  ...
)

Arguments

Details

The default is to return the fixed effects estimates together with their standard errors. If setting 'type="complex"' or to an abbreviation (e.g. "c") additional inferential statistics on these estimates will be returned, i.e. "t Value", "DF" and respective p-value "Pr > |t|". One can choose one of three denominator degrees of freedom ('ddfm')-methods. The implementation of these methods are an attempt to align with the results of SAS PROC MIXED. See the respective SAS-documentation for details.

Examples

## Not run: 
data(dataEP05A2_1)
fit <- anovaVCA(y~day/(run), dataEP05A2_1)
fixef(fit)

# for complex models it might take some time computing complex output
data(VCAdata1)
fit <- anovaMM(y~(lot+device)/(day)/(run), VCAdata1[VCAdata1$sample==2,])
fixef(fit, "c")

## End(Not run)

Description

Function calls a fast Fortran90-implementation of the SWEEP operator using the transpose of the original augmented matrix $X'X$ (see getSSQsweep). In the sweeping step, also the C matrix, needed to obtain the variance estimates from the sum of squares and the Covariance matrix of the estimates are calculated.

Usage

Fsweep(M, asgn, thresh = 1e-10, tol = 1e-10, Ncpu = 1)

Arguments

Details

This is an utility-function not intended to be called directly.

Value

(list) with eight elements:

Author(s)

Florian Dufey [email protected]

References

Goodnight, J.H. (1979), A Tutorial on the SWEEP Operator, The American Statistician, 33:3, 149-158

Description

This utility function acutally calls function 'VCAinference' first and then extracts the requested confidence interval (CI) information from the resulting object. You can specify single variance components (VC) or multiple. Not specifying any specific VC will return all.

Usage

getCI(
  obj,
  vc = NULL,
  type = c("vc", "sd", "cv"),
  tail = c("one-sided", "two-sided"),
  conf.level = 0.95,
  quiet = FALSE
)

Arguments

Author(s)

Andre Schuetzenmeister [email protected]

Examples

data(dataEP05A2_3)
fit <- remlVCA(y~day/run, dataEP05A2_3)
getCI(fit)				# will return one-sided CI for all VC
getCI(fit, type="cv")	# now on CV-scale
getCI(fit, type="cv", conf.level=.9)
# multiple row at once
getCI(fit, vc=1:3, type="cv")
getCI(fit, vc=c("total", "error"), type="cv")

Description

There are three methods implemented, which are all available in SAS PROC MIXED, namely "contain", "residual", and "satterthwaite" approximations. See the documentation of SAS PROC MIXED for details on this topic.

Usage

getDDFM(
  obj,
  L,
  ddfm = c("contain", "residual", "satterthwaite"),
  tol = 1e-12,
  method.grad = "simple",
  opt = TRUE,
  items = NULL
)

Arguments

Details

The implementation of the Satterthwaite approximation was inspired by the code of function 'calcSatterth' of R-package 'lmerTest'.

Value

(numeric) vector with the specified type of degrees of freedom

Author(s)

Andre Schuetzenmeister [email protected]

See Also

test.fixef

Description

Determine degrees of freedom for custom linear hypotheses of fixed effects or LS Means using one of three possible approximation methods.

Usage

getDF(obj, L, method = c("contain", "residual", "satterthwaite"), ...)

Arguments

Details

This is a convenience function to determine the DFs for linear hypotheses in the same way as function test.fixef. Only the "DF" part is returned here which can be passed to other functions expecting DFs as input.

Value

(numeric) vector with the DFs for each row of 'L'

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 
data(VCAdata1)
tmpDat <- VCAdata1[VCAdata1$sample==1,]
tmpDat <- tmpDat[-c(11,51,73:76),]
fitMM <- anovaMM(y~(lot+device)/(day)/(run), tmpDat)
fitMM
L <- getL(fitMM, c("lot1-lot2", "device1-device2"))
getDF(fitMM, L)						# method="contain" is Default
getDF(fitMM, L, method="res")

getDF(fitMM, L, method="satt")		# takes quite long for this model

## End(Not run)

Description

Compute variance covariance matrix of variance components of a linear mixed model via the method stated in Giesbrecht and Burns (1985).

Usage

getGB(obj, tol = 1e-12)

Arguments

Details

This function is not intended to be called by users and therefore not exported.

Value

(matrix) corresponding to the Giesbrecht & Burns approximation of the variance-covariance matrix of variance components

Author(s)

Andre Schuetzenmeister [email protected], Florian Dufey [email protected]

References

Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York

Giesbrecht, F.G. and Burns, J.C. (1985), Two-Stage Analysis Based on a Mixed Model: Large-Sample Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486

See Also

vcovVC, remlVCA, remlMM

Examples

## Not run: 
data(dataEP05A2_3)
fit <- anovaVCA(y~day/run, dataEP05A2_3)
fit <- solveMME(fit)		# some additional matrices required
getGB(fit)

## End(Not run)

Description

Intermediate precision in this context here means any sum of variances below the full model originally fitted. A typical use case could be reproducibility-experiments with a single lot or multiple lots, where a pooled version of within-lab precision shall be determined.

Usage

getIP.remlVCA(obj, vc)

Arguments

Author(s)

Andre Schuetzenmeister [email protected]

Examples

data(dataEP05A2_3)
res <- remlVCA(y~day/run, dataEP05A2_3)
IPday <- getIP.remlVCA(res, "day:run")
VCAinference(IPday)

Description

Function constructs coefficient/contrast matrices from a string-representation of linear hypotheses.

Usage

getL(obj, s, what = c("fixef", "lsmeans"))

Arguments

Details

Function constructs matrices expressing custom linear hypotheses of fixed effects or LS Means. The user has to specify a string denoting this contrast which is then transformed into a coefficient/contrast matrix. This string may contain names of fixed effects belonging to same same fixed term, numeric coefficients and mathematical operators "+" and "-" (see examples).

Value

(matrix) representing one linear hypothesis of fixed effects or LS Means per row

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 
data(dataEP05A2_2)
fit <- anovaMM(y~day/(run), dataEP05A2_2)
L <- getL(fit, c("day1-day2", "day5-day10"), what="fixef")
L
test.fixef(fit, L=L)

# another custom hypothesis
L2 <- getL(fit, "0.25*day1+0.25*day2+0.5*day3-0.5*day4-0.5*day5")
L2

# more complex model
data(VCAdata1)
dataS2 <- VCAdata1[VCAdata1$sample==2,]
fit.S2 <- anovaMM(y~(lot+device)/day/(run), dataS2)
L3 <- getL(fit.S2, c("lot1-lot2", "lot1:device3:day19-lot1:device3:day20", 
"lot1:device1:day1-lot1:device1:day2"))
L3
test.fixef(fit.S2, L3)

## End(Not run)

Description

For convenience only, extracting a specific matrix from the "Matrices" element of a 'VCA' object if this matrix exists.

Usage

getMat(obj, mat)

Arguments

Details

When 'mat="Z"' the design matrix of random effects will be returned. If one is interested in the design matrix of random effects for a specific variance component use a name like "Z" + NAME, where NAME has to be equal to the name of the VC in the 'VCA' object (see examples). The same applies to the A-matrices in the quadratic forms, use "A" + NAME for extracting a specific A-matrix.

Value

(matrix) as requested by the user

Examples

## Not run: 
data(dataEP05A2_1)
fit <- anovaVCA(y~day/run, dataEP05A2_1)
getMat(fit, "Z")
getMat(fit, "Zday")
getMat(fit, "Zday:run")
getMat(fit, "Zerror")
fit2 <- anovaMM(y~day/(run), dataEP05A2_1)
getMat(fit2, "V")			 	# Var(y)
getMat(fit2, "G")				# Var(re)

## End(Not run)

Description

Function getMM constructs overparameterized design matrices from a model formula and a data.frame.

Usage

getMM(form, Data, keep.order = TRUE)

Arguments

Details

This function constructs the overparameterized design matrix for a given dataset 'Data' according to the model formula 'form'. Each combination of factor-levels and or numeric variables is identified and accounted for by a separate column. See examples for differences compared to function 'model.matrix' (stats). This type of design matrix is used e.g. in constructing A-matrices of quadratic forms in $y$ expressing ANOVA sums of squares as such. This is key functionality of functions anovaVCA and anovaMM used e.g. in constructing the coefficient matrix $C$ whose inverse is used in solving for ANOVA Type-1 based variance components..

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 
# load example data (CLSI EP05-A2 Within-Lab Precision Experiment)
data(dataEP05A2_3)
tmpData <- dataEP05A2_3[1:10,] 

# check out the differences
getMM(~day+day:run, tmpData)
model.matrix(~day+day:run, tmpData)

# adapt factor variables in 'tmpData'
tmpData$day <- factor(tmpData$day)

# check out the differences now
getMM(~day+day:run, tmpData)
model.matrix(~day+day:run, tmpData)

# numeric covariate 'cov'
tmpData2 <- dataEP05A2_3[1:10,] 
tmpData2$cov <- 10+rnorm(10,,3)
model.matrix(~day*cov, tmpData2)

## End(Not run)

Description

Compute ANOVA Type-1 sum of squares for linear models.

Usage

getSSQsweep(Data, tobj, random = NULL)

Arguments

Details

This function performs estimation of ANOVA Type-1 sum of squares using the SWEEP-operator (see reference), operating on the augmented matrix $X'X$, where $X$ represents the design matrix not differentiating between fixed and random factors. See the numerical example in Fsweep exemplifying the type of augmentation of $X'X$ on which sweeping is carried out.

This is an utility function not intended to be called directly. For each term in the formula the design-matrix $Z$ is constructed. Matrix $X$ corresponds to binding all these $Z$-matrices together column-wise.

Degrees of freedom for each term are determined by subtracting the number of linearly dependent columns from the total number of column in X asigned to a specific term.

Value

(list) representing the with variables:

Author(s)

Andre Schuetzenmeister [email protected], Florian Dufey [email protected]

References

Goodnight, J.H., (1979), A Tutorial on the SWEEP Operator, The American Statistician, 33:3, p.149-158

See Also

Fsweep

Examples

## Not run: 
data(dataEP05A2_1)
res <- VCA:::getSSQsweep(dataEP05A2_1, terms(y~day/run))
str(res)

## End(Not run)

Description

Determine the estimated variance-covariance matrix of observations $y$.

Usage

getV(obj)

Arguments

Details

A linear mixed model can be written as $y = Xb + Zg + e$, where $y$ is the column vector of observations, $X$ and $Z$ are design matrices assigning fixed ($b$), respectively, random ($g$) effects to observations, and $e$ is the column vector of residual errors. The variance-covariance matrix of $y$ is equal to $Var(y) = ZGZ^{-T} + R$, where $R$ is the variance-covariance matrix of $e$ and $G$ is the variance-covariance matrix of $g$. Here, $G$ is assumed to be a diagonal matrix, i.e. all random effects $g$ are mutually independent (uncorrelated).

Value

(VCA) object with additional elements in the 'Matrices' element, including matrix $V$.

Author(s)

Andre Schuetzenmeister [email protected]

Description

This data set consists of the Glucose intermediate precision data in the CLSI EP05-A3 guideline, i.e. total variance for a fully-nested design with 3 variance components (day, run, error).
Note, that the results in the original EP05-A3 guideline were obtained using rounded intermediate results, whereas, package VCA uses full precision. Any differences between results listed in the CLSI EP05-A3 guideline and those generated by the package are due to error propagation in the working example presented in the CLSI guideline. Here, full precision is used for all intermediate results.

Usage

data(Glucose)

Format

data.frame with 80 rows and 3 variables.

References

CLSI EP05-A3 - Evaluation of Precision of Quantitative Measurement Procedures; Approved Guideline - Third Edition. CLSI

Description

This dataset was added to have a very large dataset available for the unit-test suite. It is an unbalanced dataset with three variables and 8070 observations.

Usage

data(HugeData)

Format

A data frame with 8070 observations on the following 3 variables.

• y
  
  The response variable.
  
  

• VC1
  
  Variable with 8 levels corresponding to top-level variance component.
  
  

• VC2
  
  Variable with 3920 levels corresponding to 2nd-level variance component.
  

Description

Assess whether an experimental design is balanced or not.

Usage

isBalanced(form, Data, na.rm = TRUE)

Arguments

Details

This function is for internal use only. Thus, it is not exported.

The approach taken here is to check whether each cell defined by one level of a factor are all equal or not. Here, data is either balanced or unbalanced, there is no concept of "planned unbalancedness" as discussed e.g. in Searle et al. (1992) p.4. The expanded (simplified) formula is divided into main factors and nested factors, where the latter are interaction terms. The $N$-dimensional contingency table, $N$ being the number of main factors, is checked for all cells containing the same number. If there are differences, the dataset is classified as "unbalanced". All interaction terms are tested individually. Firstly, a single factor is generated from combining factor levels of the first $(n-1)$ variables in the interaction term. The last variable occuring in the interaction term is then recoded as factor-object with $M$ levels. $M$ is the number of factor levels within each factor level defined by the first $(n-1)$ variables in the interaction term. This is done to account for the independence within sub-classes emerging from the combination of the first $(n-1)$ variables.

Value

(logical) TRUE if data is balanced, FALSE if data is unbalanced (according to the definition of balance used)

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 
data1 <- data.frame(site=gl(3,8), lot=factor(rep(c(2,3,1,2,3,1), 
rep(4,6))), day=rep(1:12, rep(2,12)), y=rnorm(24,25,1))

# not all combinations of 'site' and 'lot' in 'data1'

VCA:::isBalanced(y~site+lot+site:lot:day, data1)

# balanced design for this model

VCA:::isBalanced(y~lot+lot:day, data1)

# gets unbalanced if observation is NA

data1[1,"y"] <- NA
VCA:::isBalanced(y~lot+lot:day, data1)
VCA:::isBalanced(y~lot+lot:day, data1, FALSE)

## End(Not run)

Description

This function accepts all parameters applicable in and forwards them to function legend. There will be only made some modifications to the X-coordinate ensuring that the legend is plotted in the right margin of the graphic device. Make sure that you have reserved sufficient space in the right margin, e.g. 'plot.VFP(....., mar=c(4,5,4,10))'.

Usage

legend.m(
  x = c("center", "bottomright", "bottom", "bottomleft", "left", "topleft", "top",
    "topright", "right"),
  y = NULL,
  margin = c("right", "bottomright", "bottom", "bottomleft", "left", "topleft", "top",
    "topright"),
  offset = 0.05,
  ...
)

Arguments

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 

par( mar=c(10,10,10,10) )
plot(1, type="n", axes=FALSE, xlab="", ylab="")
box()
# add legend to different regions within the 'margin'
legend.m(margin="topleft", 		fill="black",	legend=c("topleft"))
legend.m(margin="top", 			fill="red", 	legend=c("top"))
legend.m(margin="topright", 		fill="blue",	legend=c("topright"))
legend.m(margin="right", 		fill="green",	legend=c("right"))
legend.m(margin="bottomright", 	fill="yellow",	legend=c("bottomright"))
legend.m(margin="bottom", 		fill="orange",	legend=c("bottom"))
legend.m(margin="bottomleft", 	fill="cyan",	legend=c("bottomleft"))
legend.m(margin="left", 			fill="magenta", legend=c("left"))

data(dataEP05A2_3)
dataEP05A2_3$user <- sample(rep(c(1,2), 40))
 
varPlot( y~day+day:run, dataEP05A2_3, mar=c(1,5,1,7), VCnam=list(side=4), 
	        Points=list(pch=list(var="user", pch=c(2, 8))) )
# always check order of factor levels before annotating
order(unique(dataEP05A2_3$user))
legend.m(pch=c(8,2), legend=c("User 1", "User 2"))

# using different colors 
varPlot( y~day+day:run, dataEP05A2_3, mar=c(1,5,1,7), VCnam=list(side=4),
         Points=list(col=list(var="user", col=c("red", "green"))) )
legend.m(fill=c("green", "red"), legend=c("User 1", "User 2"))

# two additional classification variables
dataEP05A2_3$user <- sample(rep(c(1,2), 40))
dataEP05A2_3$cls2 <- sample(rep(c(1,2), 40))

# now combine point-coloring and plotting symbols
# to indicate two additional classification variables
varPlot( y~day+day:run, dataEP05A2_3, mar=c(1,5,1,10),
         VCnam=list(side=4, cex=1.5),
         Points=list(col=list(var="user", col=c("red", "darkgreen")),
                     pch=list(var="cls2", pch=c(21, 22)),
                     bg =list(var="user", bg =c("orange", "green"))) )

# add legend to (right) margin
legend.m(margin="right", pch=c(21, 22, 22, 22), 
         pt.bg=c("white", "white", "orange", "green"),  
         col=c("black", "black", "white", "white"), 
         pt.cex=c(1.75, 1.75, 2, 2),
         legend=c("Cls2=1", "Cls2=2", "User=2", "User=1"), 
         cex=1.5)

## End(Not run)

Description

This function restricts the variance-covariance matrix of random effects $G$ to be either diagonal ('cov=FALSE') or to take any non-zero covariances into account (default, 'cov=TRUE').

Usage

lmerG(obj, cov = FALSE)

Arguments

Details

This function is not intended to be called directly by users and therefore not exported!

Value

(Matrix) representing the variance-covariance structure of random effects $G$

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 
library(lme4)
data(Orthodont)
Ortho <- Orthodont
Ortho$age2 <- Ortho$age - 11
Ortho$Subject <- factor(as.character(Ortho$Subject))
fit <-lmer(distance~Sex+Sex:age2+(age2|Subject), Ortho) 
G1 <- VCA:::lmerG(fit, cov=FALSE)
G2 <- VCA:::lmerG(fit, cov=TRUE)
G1[1:10,1:10]
G2[1:10,1:10]

## End(Not run)

Description

Function derives and computes all matrices required for down-stream analyses of VCA-objects fitted with REML via function lmer.

Usage

lmerMatrices(obj, tab = NULL, terms = NULL, cov = FALSE, X = NULL)

Arguments

Details

Mixed Model Equations (MME) are solved for fixed and random effects applying the same constraints as in anovaMM. The most elaborate and therefore time consuming part is to prepare all matrices required for approximating the variance-covariance matrix of variance components (see getGB). To reduce the computational time, this function tries to optimize object-classes depending on whether Intel's (M)ath (K)ernel (L)ibrary could be loaded or not. MKL appears to be more performant with ordinary matrix-objects, whereas all other computations are perfomred using matrix-representations of the Matrix-package.

This function is not intended to be called directly by users and therefore not exported.

Value

(list), a premature 'VCA' object

Author(s)

Andre Schuetzenmeister [email protected]

See Also

remlVCA, remlMM

Derive VCA-Summary Table from an Object Fitted via Function lmer

Description

This function builds a variance components analysis (VCA) table from an object representing a model fitted by function lmer of the lme4 R-package.

Usage

lmerSummary(
  obj,
  VarVC = TRUE,
  terms = NULL,
  Mean = NULL,
  cov = FALSE,
  X = NULL,
  tab.only = FALSE
)

Arguments

Details

It applies the approximation of the variance-covariance matrix of variance components according to Giesbrecht & Burns (1985) and uses this information to approximate the degrees of freedom according to Satterthwaite (see SAS PROC MIXED documentation option 'CL').

This function can be used to create a VCA-results table from almost any fitted 'lmerMod'-object, i.e. one can apply it to a model fitted via function lmer of the lme4-package. The only additional argument that needs to be used is 'tab.only' (see examples).

Value

(list) still a premature 'VCA'-object but close to a "complete" 'VCA'-object

Author(s)

Andre Schuetzenmeister [email protected]

References

Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York

Giesbrecht, F.G. and Burns, J.C. (1985), Two-Stage Analysis Based on a Mixed Model: Large-Sample Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486

See Also

remlVCA, remlMM

Examples

## Not run: 
# fit a model with a VCA-function first
data(VCAdata1)
fit0 <- remlVCA(y~(device+lot)/day/run, subset(VCAdata1, sample==5))
fit0

# fit the same model with function 'lmer' of the 'lme4'-package
library(lme4)
fit1 <- lmer(y~(1|device)+(1|lot)+(1|device:lot:day)+(1|device:lot:day:run),
subset(VCAdata1, sample==5))
lmerSummary(fit1, tab.only=TRUE)

## End(Not run)

Description

This function is taken from the Rprofile.site file of Microsoft R Open. It was added to the package namespace to avoid a NOTE during the R CMD check process stating that this function is not gobally defined.

Usage

load_if_installed(package)

Arguments

Details

Only change to the original version is a different bracketing scheme to match the one used in the remaining source-code of the package.

Author(s)

Authors of the Rprofile.site file in Microsoft R Open.

Description

Computes Least Squares Means (LS Means) of fixed effects for fitted mixed models of class 'VCA'.

Usage

lsmeans(
  obj,
  var = NULL,
  type = c("simple", "complex"),
  ddfm = c("contain", "residual", "satterthwaite"),
  at = NULL,
  contr.mat = FALSE,
  quiet = FALSE
)

Arguments

Details

Function computes LS Means of fixed effects and their corresponding standard errors. In case of setting argument 'type' equal to "complex" (or any abbreviation) a $t$-test is performed on each LS Mean, returning degrees of freedom, t-statistic and corresponding p-values. One can choose from one of three denominator degrees of freedom ('ddfm')-methods.

Actually, function test.fixef is called with the "no intercept" version of the fitted model. The "complex" option is significantly slower for unbalanced designs (see test.fixef for details). In case that the 'VarCov' element of the 'VCA' object already exists (calling vcovVC), which is the most time consuming part, results can be obtained in less amount of time.

Standard Errors of LS Means are computed as $TPT^{T}$, where $T$ is the LS Means generating contrast matrix and $P$ is the variance-covariance matrix of fixed effects.

Argument at can be used to modify the values of covariables when computing LS Means and/or to apply different weighting schemes for (fixed) factor variables in the model, e.g. when the prevelance of factor-levels differs from a uniform distribution. Usually, if the weighting scheme is not modified, each factor-level will contribute $1/N$ to the LS Mean, where $N$ corresponds to the number of factor-levels.

Covariables have to be specified as 'name=value', where value can be a vector of length > 1. Each value will be evaluated for each row of the original LS Means contrast matrix. If multiple covariables are specified, the i-th element of covariable 1 will be matched with the i-th element of covariable(s) 2...M, where $M$ is the number of covariables in the model.

To apply a different weighting scheme for factor-variables one has to specify 'factor-name=c(level-name_1=value_1, level-name_2=value_2, ..., level-name_N=value_N)'. The sum of all 'value_i' elements must be equal to 1, otherwise, this factor-variable will be skipped issuing a warning. If any levels 'level-name_i' cannot be found for factor-variable 'factor-name', this variable will also be skipped and a warning will be issued. See the examples section to get an impression of how this works.

Value

(matrix) with LS Means of fixed effects and respective standard errors, in case of 'type="complex"'

Author(s)

Andre Schuetzenmeister [email protected]

Examples

#
## Not run: 
data(dataEP05A2_2)
fit1 <- anovaMM(y~day/(run), dataEP05A2_2)
lsmeans(fit1)
lsmeans(fit1,, "complex")

# a more complex model
data(VCAdata1)
fit2 <- anovaMM(y~(lot+device)/(day)/(run), VCAdata1[VCAdata1$sample==2,])
lsmeans(fit2, "lot")
lsmeans(fit2, "device", "complex")

# pre-computed 'VarCov' element saves time
system.time(lsm1 <- lsmeans(fit2, "device", "complex"))
fit2$VarCov <- vcovVC(fit2)
system.time(lsm2 <- lsmeans(fit2, "device", "complex"))
lsm1
lsm2 

# simulate some random data 
set.seed(212)
id <- rep(1:10,10)
x <- rnorm(200)
time <- sample(1:5,200,replace=T)
y <- rnorm(200)+time
snp <- sample(0:1,200,replace=T)
dat <- data.frame(id=id,x=x,y=y,time=time,snp=snp)
dat$snp <- as.factor(dat$snp)
dat$id <- as.factor(dat$id)
dat$time <- as.numeric(dat$time)
dat$sex <- gl(2, 100, labels=c("Male", "Female"))
dat$y <- dat$y + rep(rnorm(2, 5, 1), c(100, 100))

fit3 <- remlMM(y~snp+time+snp:time+sex+(id)+(id):time, dat)

# compute standard LS Means for variable "snp"
lsmeans(fit3, var="snp")
lsmeans(fit3, var="snp", type="c")    # comprehensive output

# compute LS Means at timepoints 1, 2, 3, 4
# Note: original LS Means are always part of the output
lsmeans(fit3, var="snp", at=list(time=1:4))

# compute LS Means with different weighting scheme
# for factor-variable 'sex'
lsmeans(fit3, var="snp", at=list(sex=c(Male=.3, Female=.7)))

# combine covariables at some value and altering the
# weighting scheme
lsmeans(fit3, var="snp", at=list(time=1:4, sex=c(Male=.3, Female=.7)))

# now with comprehensive output and requesting the
# LS Means generating contrast matrix
lsmeans(fit3, var="snp", type="complex", contr.mat=TRUE,
at=list(time=1:4, sex=c(Male=.3, Female=.7)))

## End(Not run)

Description

This data set is used for unit-testing LS Means functionality. Reference results were generated in SAS PROC MIXED and rounded to two decimals as covariance parameter estimates slightly differ.

Usage

data(LSMeans_Data)

Format

data.frame with 200 rows and 6 variables.

Description

Function determines appropriate contrast matrix for computing the LS Means of each factor level of one or multiple fixed effects variables.

Usage

lsmMat(obj, var = NULL, quiet = FALSE)

Arguments

Details

This functions implements the 5 rules given in the documentation of SAS PROC GLM for computing the LS Means.#' The LS Means correspond to marginal means adjusted for bias introduced by unbalancedness.

Value

(matrix) where each row corresponds to a LS Means generating contrast for each factor level of one or multiple fixed effects variable(s)

Author(s)

Andre Schutzenmeister [email protected]

Examples

## Not run: 
data(dataEP05A2_1)
fit1 <- anovaMM(y~day/run, dataEP05A2_1)

VCA:::lsmMat(fit1, "day")	# function not exported
VCA:::lsmMat(fit1, "run")
VCA:::lsmMat(fit1)			# is equal to listing all fixed terms

# a more complex and unbalanced model
data(VCAdata1)
datS1 <- VCAdata1[VCAdata1$sample == 1, ]
set.seed(42)
datS1ub <- datS1[-sample(1:nrow(datS1))[1:25],]
fit2 <- anovaMM(y~(lot+device)/day/(run), datS1ub)
VCA:::lsmMat(fit2, c("lot", "device"))

## End(Not run)

Description

This data set consists of 754 observations. There are 3 laboratories (Lab), 3 lots (Lot), 21 days (Days) per lab-lot combination, and 2 runs per day. The response variable is Result. This dataset is used in examples and unit-tests (see subdir 'UnitTests' of the package-dir).

Usage

data(MLrepro)

Format

data.frame with 754 rows and 5 variables.

Description

Function returns the data-element of 'object' and adds the terms-element as attribute.

Usage

## S3 method for class 'VCA'
model.frame(formula, ...)

Arguments

Details

It enables application of functions relying on the existence of this method, e.g. the functin 'glht' of the 'multcomp' R-package.

Value

(data.frame) with attribute 'terms'

Author(s)

Andre Schuetzenmeister [email protected]

Description

Function returns matrix X corresponding to the design matrix of fixed effects of the fitted model.

Usage

## S3 method for class 'VCA'
model.matrix(object, ...)

Arguments

Description

This function is originally impelemented in package 'MASS' as function ginv. It was adapted to be able to deal with matrices from the 'Matrix' package, e.g. sparse matrices.

Usage

MPinv(X, tol = sqrt(.Machine$double.eps))

Arguments

Value

(object) A Moore-Penrose inverse of X.

Author(s)

Authors of the 'MASS' package.

Description

Functions attempts to standardize input data for linear mixed model analyses to overcome the problem that analysis results sometimes depend on ordering of the data and definition of factor-levels.

Usage

orderData(Data, trms, order.data = TRUE, exclude.numeric = TRUE, quiet = FALSE)

Arguments

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 
# random ordering
data(dataEP05A2_1)
dat <- dataEP05A2_1
levels(dat$day) <- sample(levels(dat$day))
# this has direct impact e.g. on order of estimated effects
fit <- anovaVCA(y~day/run, dat, order.data=FALSE)
ranef(fit)
# to guarantee consistent analysis results
# independent of the any data orderings option
# 'order.data' is per default set to TRUE:
fit <- anovaVCA(y~day/run, dat)
ranef(fit)
# which is identical to:
fit2 <- anovaVCA(y~day/run, orderData(dat, y~day/run), order.data=FALSE)
ranef(fit2)

## End(Not run)

Description

The Orthodont data frame has 108 rows and 4 columns of the change in an orthdontic measurement over time for several young subjects.

This dataset was included to simplify its usage in automated unit-tests (see directory UnitTests) and examples.

Investigators at the University of North Carolina Dental School followed the growth of 27 children (16 males, 11 females) from age 8 until age 14. Every two years they measured the distance between the pituitary and the pterygomaxillary fissure, two points that are easily identified on x-ray exposures of the side of the head.

Usage

data(Orthodont)

Format

data.frame with 80 rows and 3 variables.

References

Pinheiro, J. C. and Bates, D. M. (2000), Mixed-Effects Models in S and S-PLUS, Springer, New York. (Appendix A.17)

Potthoff, R. F. and Roy, S. N. (1964), A generalized multivariate analysis of variance model useful especially for growth curve problems, Biometrika, 51, 313-326.

Description

Create a variability chart from a 'VCA'-object, i.e. from a fitted model.

Usage

## S3 method for class 'VCA'
plot(x, ...)

Arguments

Details

This function extracts the data and the model-formula from a fitted 'VCA'-object and calls function varPlot accepting all its arguments. Please see the documention of function varPlot for a detailed description.

It will not be differentiated between fixed and random effects when calling this function on a fitted linear mixed model.

Value

nothing, instead a plot is generated

Author(s)

Andre Schuetzenmeister [email protected]

See Also

varPlot, anovaVCA,remlVCA, anovaMM,remlMM

Examples

## Not run: 
data(dataEP05A2_1)
fit <- anovaVCA(y~day/run, dataEP05A2_1)

# standard plot without any extras
plot(fit)

# plot with some additional features
plot(fit, MeanLine=list(var=c("int", "day"), col=c("cyan", "blue"), lwd=c(2,2)))

# more complex model
data(realData)
Data <- realData[realData$PID == 1,]
fit2 <- anovaVCA(y~(calibration+lot)/day/run, Data)
plot(fit2, 
		BG=list(var="calibration",
				col=c("#f7fcfd","#e5f5f9","#ccece6","#99d8c9",
				      "#66c2a4","#41ae76","#238b45","#006d2c","#00441b"),
				col.table=TRUE),
	VLine=list(var=c("calibration", "lot"),
				   col=c("black", "darkgray"), lwd=c(2,1), col.table=TRUE),
		JoinLevels=list(var="lot", col=c("#ffffb2","orangered","#feb24c"),
				        lwd=c(2,2,2)),
		MeanLine=list(var="lot", col="blue", lwd=2))


## End(Not run)

Description

Plots, possibly transformed, random variates of a linear mixed model (random effects, contitional or marginal residuals).

Usage

plotRandVar(
  obj,
  term = NULL,
  mode = c("raw", "student", "standard", "pearson"),
  main = NULL,
  Xlabels = list(),
  Points = list(),
  Vlines = list(),
  pick = FALSE,
  ...
)

Arguments

Details

Function plots either random effects of a 'VCA' object or residuals. Parameter 'term' is used to specify either one. If 'term' is one of c("conditional", "marginal") corresponding residuals will be plotted (see resid for details). If 'term' is either the name of a random term in the formula of the 'VCA' object or an integer specifying the i-th random term, corresponding random effects will be plotted. Both types of random variates (random effects, residuals) can be plotted untransformed ("raw"), "studentized" or "standardized". In case of residuals, one can also use the "Pearson"-type transformation.

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 
data(dataEP05A2_1)
fit <- anovaVCA(y~day/run, dataEP05A2_1)
# solve mixed model equations including random effects
fit <- solveMME(fit)
plotRandVar(fit, "cond", "stand")	
plotRandVar(fit, 1, "stud")						# 1st random term 'day'
plotRandVar(fit, "day", "stud")					# equivalent to the above

# for larger datasets residuals can hardly be identified
# pick out interesting points with the mouse

plotRandVar(fit, "marg", "stud", pick=TRUE)

# customize the appearance
plotRandVar( fit, 1, "stud", Vlines=list(col=c("red", "darkgreen")), 
	Xlabels=list(offset=.5, srt=60, cex=1, col="blue"),
	Points=list(col=c("black", "red", rep("black", 18)),
pch=c(3,17,rep(3,18)), cex=c(1,2,rep(1,18))))	

## End(Not run)

Description

Model returns fitted values in case newdata is NULL or evaluates the fitted model employing user-specified data newdata. The default is that fitted values incorporate fixed effects and random effects, leaving out the (conditional) residuals only. If the interest lies in constraining predictions to the fixed effects only set re=NA or incorporate just part of the random variability specifying distinct random effects (see re.

Usage

## S3 method for class 'VCA'
predict(object, newdata = NULL, re = NULL, allow.new.levels = FALSE, ...)

Arguments

Value

(numeric) vector of prediction results

Author(s)

Andre Schuetzenmeister [email protected]

Examples

## Not run: 
# fit LMM with fixed lot and device effects and test for lot-differences
data(VCAdata1)
datS5 <- subset(VCAdata1, sample==5)
fitS5 <- fitLMM(y~(lot+device)/(day)/(run), datS5, "anova")
fitS5

# fitted values including fixed and random effects
pred0 <- predict(fitS5)
pred0
# sanity check:
all(round(pred0 + resid(fitS5) - datS5$y, 12) == 0)
# restrict to fixed effects
predict(fitS5, re=NA)
# restrict to fixed effects and dayly random effects
# see required names
fitS5$random
predict(fitS5, re="lot:device:day")

# check against original 'lmer'-predictions
# use version from VCA-package (ordinary data.frame)
data(Orthodont, package="VCA")
Ortho <- Orthodont
Ortho$age2 <- Ortho$age-11
# use exactly the same data, same ordering
Ortho <- orderData(Ortho, distance ~ Sex * age2 + (Subject) * age2)
fit.fitLMM <- fitLMM(distance ~ Sex * age2 + (Subject) * age2, Ortho, "reml")
library(lme4)
fit.lmer <- lmer(distance ~ Sex + age2 + Sex:age2 + (age2 | Subject), Ortho)
# check fitted value first (fixed + random effects)
predict(fit.lmer)
predict(fit.fitLMM)
# restrict to fixed part only
predict(fit.lmer, re.form=NA)  
predict(fit.fitLMM, re=NA)
# user-specified 'newdata'
newdata <- Ortho[45:54,]
newdata$age2 <- newdata$age2 + 5
# include fixed and random effects
predict(fit.lmer, newdata)
predict(fit.fitLMM, newdata)
# generate new data
newdata <- Ortho[45:54,]          
newdata$age2 <- newdata$age2 + 5
# predict on newdata using fixed and random effects
predict(fit.lmer, newdata) 
predict(fit.fitLMM, newdata)       
# restrict prediction to random Subject effects
predict(fit.lmer, newdata, re.form=~(1|Subject))        
predict(fit.fitLMM, newdata, re="Subject")
# restrict prediction to random per-Subject slope
predict(fit.lmer, newdata, re.form=~(age2-1|Subject)) 
predict(fit.fitLMM, newdata, re="age2:Subject")

## End(Not run)