This is a first attempt at an R version of a simple nested model. This is the beginning of the R code/analyses comparable to that for the SAS nested/sub-sampling example, using 3 treatments and 4 trees per treatment, and weighing 6 apples per tree.
library(lmerTest)
library(doBy)
# read data, fit Y = mu + trt + tree(trt) + e
# gives same results as SAS proc mixed (good)
# set a string as the path/filename, so we can easily move the
# file from one directory to another
# this will typically be (for a Windows system, something
# like "c:/Users/UserName/wherever_your_R_data_are/"
# for a Linux system it will probably be "/home/UserName/your_R_files/'
FilePath <- "complete path to by data subdirectory/"
FileName <- "datafile_name"
filen <- paste(FilePath,FileName,sep="")
filen
ds<-read.table(filen,header=TRUE)
ds$tree1 <- (ds$trt-1)*4 + ds$tree
ds$trt<-as.factor(ds$trt)
ds$tree1 <- as.factor(ds$tree1)
# get basic summary of our dataframe
summary(ds)
# fit mixed model, trt and tree nested within trt
lm1 <- lmer(wt ~ trt + (1|trt:tree1),data=ds)
str(lm1)
lm1
# output fixed effects estimate/solution vector
# NB just estimates, no se
fixef(lm1)
anova(lm1)
summary(lm1)
rand(lm1)
# lmerTest has methods for lsmeans and pairwise differences
# we can use esticon for other, roll-your-own tests
#
lsmeans(lm1)
difflsmeans(lm1)
# using the esticon function from package doBy,
# we can explicitly calculate our own lsmeans or whatever
# this can be useful if we want/need to calculate our
# own estimates of various things, and then use Scheffe's test
w1 <- c(1,0,0)
esticon(lm1,w1)
w1 <- c(1,1,0)
esticon(lm1,w1)
w1 <- c(1,0,1)
esticon(lm1,w1)
# consider that we need to test the average of
# treatments 1 and 2 vs trt 3. This is not provided
# by either the lsmeans or the difflsmeans.
# construct our own estimate/contrast
w1 <- c(0,0.5,-1)
diff1 <- esticon(lm1,w1)
diff1
# get tabulated F value, df1 = ndf = Number of treatments - 1, = (3-1=2)
# df2 = ddf = N trees - N treatments (12-3 = 9)
Ftab <- qf(0.95,df1=2,df2=9)
Ftab
ndf <- 2
CD <- diff1$Std.Error * sqrt(ndf*Ftab)
CD
# Estimate = -34.23, se = 4.8449
# CD = 14.1361
# |Estimate| > CD, therefore SS
#
Results, raw
> library(lmerTest)
Loading required package: Matrix
Loading required package: lme4
Attaching package: lmerTest
The following object is masked from package:lme4:
lmer
The following object is masked from package:stats:
step
>
> library(doBy)
Loading required package: survival
>
> # read data, fit Y = mu + trt + tree(trt) + e
> # gives same results as SAS proc mixed (good)
>
> # set a string as the path/filename, so we can easily move the
> # file from one directory to another
> # this will typically be (for a Windows system, something
> # like "c:/Users/UserName/wherever_your_R_data_are/"
> # for a Linux system it will probably be "/home/UserName/your_R_files/'
>
> FilePath <- "c:/Users/Roger Cue/wwwstats/R/nested/"
> FileName <- "appletree1.txt"
>
> filen <- paste(FilePath,FileName,sep="")
> filen
[1] "c:/Users/Roger Cue/wwwstats/R/nested/appletree1.txt"
>
> ds<-read.table(filen,header=TRUE)
> ds$tree1 <- (ds$trt-1)*4 + ds$tree
> ds$trt<-as.factor(ds$trt)
> ds$tree1 <- as.factor(ds$tree1)
>
> # get basic summary of our dataframe
> summary(ds)
trt tree apple wt tree1
1:24 Min. :1.00 Min. :1.0 Min. :313.1 1 : 6
2:24 1st Qu.:1.75 1st Qu.:2.0 1st Qu.:334.9 2 : 6
3:24 Median :2.50 Median :3.5 Median :349.3 3 : 6
Mean :2.50 Mean :3.5 Mean :351.3 4 : 6
3rd Qu.:3.25 3rd Qu.:5.0 3rd Qu.:365.1 5 : 6
Max. :4.00 Max. :6.0 Max. :391.8 6 : 6
(Other):36
>
> # fit mixed model, trt and tree nested within trt
> lm1 <- lmer(wt ~ trt + (1|trt:tree1),data=ds)
>
> str(lm1)
Formal class 'merModLmerTest' [package "lmerTest"] with 13 slots
..@ Gp : int [1:2] 0 12
..@ call : language lme4::lmer(formula = wt ~ trt + (1 | trt:tree1), data = ds)
..@ frame :'data.frame': 72 obs. of 3 variables:
.. ..$ wt : num [1:72] 313 329 334 330 335 ...
.. ..$ trt : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ tree1: Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 1 1 2 2 2 2 ...
.. ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 wt ~ trt + (1 + trt:tree1)
.. .. .. ..- attr(*, "variables")= language list(wt, trt, tree1)
.. .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 2 1
.. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. ..$ : chr [1:3] "wt" "trt" "tree1"
.. .. .. .. .. ..$ : chr [1:2] "trt" "trt:tree1"
.. .. .. ..- attr(*, "term.labels")= chr [1:2] "trt" "trt:tree1"
.. .. .. ..- attr(*, "order")= int [1:2] 1 2
.. .. .. ..- attr(*, "intercept")= int 1
.. .. .. ..- attr(*, "response")= int 1
.. .. .. ..- attr(*, ".Environment")=
.. .. .. ..- attr(*, "predvars")= language list(wt, trt, tree1)
.. .. .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "factor" "factor"
.. .. .. .. ..- attr(*, "names")= chr [1:3] "wt" "trt" "tree1"
.. .. .. ..- attr(*, "predvars.fixed")= language list(wt, trt)
.. ..- attr(*, "formula")=Class 'formula' length 3 wt ~ trt + (1 | trt:tree1)
.. .. .. ..- attr(*, ".Environment")=
..@ flist :List of 1
.. ..$ trt:tree1: Factor w/ 12 levels "1:1","1:2","1:3",..: 1 1 1 1 1 1 2 2 2 2 ...
.. ..- attr(*, "assign")= int 1
..@ cnms :List of 1
.. ..$ trt:tree1: chr "(Intercept)"
..@ lower : num 0
..@ theta : num 0.717
..@ beta : num [1:3] 333.5 12.9 40.7
..@ u : num [1:12] -6.007 0.866 7.471 -2.33 0.444 ...
..@ devcomp:List of 2
.. ..$ cmp : Named num [1:10] 16.88 5.31 5723.57 625.23 6348.8 ...
.. .. ..- attr(*, "names")= chr [1:10] "ldL2" "ldRX2" "wrss" "ussq" ...
.. ..$ dims: Named int [1:12] 72 72 3 69 1 12 1 1 0 3 ...
.. .. ..- attr(*, "names")= chr [1:12] "N" "n" "p" "nmp" ...
..@ pp :Reference class 'merPredD' [package "lme4"] with 18 fields
.. ..$ Lambdat:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. ..@ i : int [1:12] 0 1 2 3 4 5 6 7 8 9 ...
.. .. .. ..@ p : int [1:13] 0 1 2 3 4 5 6 7 8 9 ...
.. .. .. ..@ Dim : int [1:2] 12 12
.. .. .. ..@ Dimnames:List of 2
.. .. .. .. ..$ : NULL
.. .. .. .. ..$ : NULL
.. .. .. ..@ x : num [1:12] 0.717 0.717 0.717 0.717 0.717 ...
.. .. .. ..@ factors : list()
.. ..$ LamtUt :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. ..@ i : int [1:72] 0 0 0 0 0 0 1 1 1 1 ...
.. .. .. ..@ p : int [1:73] 0 1 2 3 4 5 6 7 8 9 ...
.. .. .. ..@ Dim : int [1:2] 12 72
.. .. .. ..@ Dimnames:List of 2
.. .. .. .. ..$ : NULL
.. .. .. .. ..$ : chr [1:72] "1" "2" "3" "4" ...
.. .. .. ..@ x : num [1:72] 0.717 0.717 0.717 0.717 0.717 ...
.. .. .. ..@ factors : list()
.. ..$ Lind : int [1:12] 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ Ptr :
.. ..$ RZX : num [1:12, 1:3] 2.13 2.13 2.13 2.13 2.13 ...
.. ..$ Ut :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. ..@ i : int [1:72] 0 0 0 0 0 0 1 1 1 1 ...
.. .. .. ..@ p : int [1:73] 0 1 2 3 4 5 6 7 8 9 ...
.. .. .. ..@ Dim : int [1:2] 12 72
.. .. .. ..@ Dimnames:List of 2
.. .. .. .. ..$ : chr [1:12] "1:1" "1:2" "1:3" "1:4" ...
.. .. .. .. ..$ : chr [1:72] "1" "2" "3" "4" ...
.. .. .. ..@ x : num [1:72] 1 1 1 1 1 1 1 1 1 1 ...
.. .. .. ..@ factors : list()
.. ..$ Utr : num [1:12] 1409 1438 1464 1424 1491 ...
.. ..$ V : num [1:72, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ VtV : num [1:3, 1:3] 72 0 0 24 24 0 24 0 24
.. ..$ Vtr : num [1:3] 25296 8313 8980
.. ..$ X : num [1:72, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:72] "1" "2" "3" "4" ...
.. .. .. ..$ : chr [1:3] "(Intercept)" "trt2" "trt3"
.. .. ..- attr(*, "assign")= int [1:3] 0 1 1
.. .. ..- attr(*, "contrasts")=List of 1
.. .. .. ..$ trt: chr "contr.treatment"
.. .. ..- attr(*, "msgScaleX")= chr(0)
.. ..$ Xwts : num [1:72] 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ Zt :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. ..@ i : int [1:72] 0 0 0 0 0 0 1 1 1 1 ...
.. .. .. ..@ p : int [1:73] 0 1 2 3 4 5 6 7 8 9 ...
.. .. .. ..@ Dim : int [1:2] 12 72
.. .. .. ..@ Dimnames:List of 2
.. .. .. .. ..$ : chr [1:12] "1:1" "1:2" "1:3" "1:4" ...
.. .. .. .. ..$ : chr [1:72] "1" "2" "3" "4" ...
.. .. .. ..@ x : num [1:72] 1 1 1 1 1 1 1 1 1 1 ...
.. .. .. ..@ factors : list()
.. ..$ beta0 : num [1:3] 0 0 0
.. ..$ delb : num [1:3] 333.5 12.9 40.7
.. ..$ delu : num [1:12] -6.007 0.866 7.471 -2.33 0.444 ...
.. ..$ theta : num 0.717
.. ..$ u0 : num [1:12] 0 0 0 0 0 0 0 0 0 0 ...
.. ..and 43 methods, of which 31 are possibly relevant:
.. .. b, beta, CcNumer, copy#envRefClass, initialize, initializePtr,
.. .. installPars, L, ldL2, ldRX2, linPred, P, ptr, RX, RXdiag, RXi, setBeta0,
.. .. setDelb, setDelu, setTheta, setZt, solve, solveU, sqrL, u, unsc,
.. .. updateDecomp, updateL, updateLamtUt, updateRes, updateXwts
..@ optinfo:List of 7
.. ..$ optimizer: chr "bobyqa"
.. ..$ control :List of 1
.. .. ..$ iprint: int 0
.. ..$ derivs :List of 2
.. .. ..$ gradient: num 7.67e-08
.. .. ..$ Hessian : num [1, 1] 34.7
.. ..$ conv :List of 2
.. .. ..$ opt : int 0
.. .. ..$ lme4: list()
.. ..$ feval : int 19
.. ..$ warnings : list()
.. ..$ val : num 0.717
..@ resp :Reference class 'lmerResp' [package "lme4"] with 9 fields
.. ..$ Ptr :
.. ..$ mu : num [1:72] 329 329 329 329 329 ...
.. ..$ offset : num [1:72] 0 0 0 0 0 0 0 0 0 0 ...
.. ..$ sqrtXwt: num [1:72] 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ sqrtrwt: num [1:72] 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ weights: num [1:72] 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ wtres : num [1:72] -16.1045 -0.0355 5.1105 0.9205 5.8195 ...
.. ..$ y : num [1:72] 313 329 334 330 335 ...
.. ..$ REML : int 3
.. ..and 26 methods, of which 14 are possibly relevant:
.. .. allInfo, copy#envRefClass, initialize, initialize#lmResp, initializePtr,
.. .. initializePtr#lmResp, objective, ptr, ptr#lmResp, setOffset, setResp,
.. .. setWeights, updateMu, wrss
> lm1
Linear mixed model fit by REML ['merModLmerTest']
Formula: wt ~ trt + (1 | trt:tree1)
Data: ds
REML criterion at convergence: 530.0187
Random effects:
Groups Name Std.Dev.
trt:tree1 (Intercept) 6.875
Residual 9.592
Number of obs: 72, groups: trt:tree1, 12
Fixed Effects:
(Intercept) trt2 trt3
333.47 12.90 40.68
>
> # output fixed effects estimate/solution vector
> # NB just estimates, no se
> fixef(lm1)
(Intercept) trt2 trt3
333.47283 12.89762 40.68225
>
> anova(lm1)
Analysis of Variance Table of type III with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
trt 5082.8 2541.4 2 9 27.62 0.0001442 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
> summary(lm1)
Linear mixed model fit by REML t-tests use Satterthwaite approximations to
degrees of freedom [lmerMod]
Formula: wt ~ trt + (1 | trt:tree1)
Data: ds
REML criterion at convergence: 530
Scaled residuals:
Min 1Q Median 3Q Max
-1.90165 -0.73876 0.08265 0.62441 1.99659
Random effects:
Groups Name Variance Std.Dev.
trt:tree1 (Intercept) 47.26 6.875
Residual 92.01 9.592
Number of obs: 72, groups: trt:tree1, 12
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 333.473 3.956 9.000 84.298 2.35e-14 ***
trt2 12.898 5.594 9.000 2.305 0.0466 *
trt3 40.682 5.594 9.000 7.272 4.70e-05 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Correlation of Fixed Effects:
(Intr) trt2
trt2 -0.707
trt3 -0.707 0.500
> rand(lm1)
Analysis of Random effects Table:
Chi.sq Chi.DF p.value
trt:tree1 10.7 1 0.001 **
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
>
> # lmerTest has methods for lsmeans and pairwise differences
> # we can use esticon for other, roll-your-own tests
> #
> lsmeans(lm1)
Least Squares Means table:
trt Estimate Standard Error DF t-value Lower CI Upper CI p-value
trt 1 1 333.47 3.96 9 84.30 325 342 <2e-16 ***
trt 2 2 346.37 3.96 9 87.56 337 355 <2e-16 ***
trt 3 3 374.16 3.96 9 94.58 365 383 <2e-16 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
> difflsmeans(lm1)
Differences of LSMEANS:
Estimate Standard Error DF t-value Lower CI Upper CI p-value
trt 1 - 2 -12.9 5.59 9.0 -2.31 -25.6 -0.242 0.05 *
trt 1 - 3 -40.7 5.59 9.0 -7.27 -53.3 -28.027 <2e-16 ***
trt 2 - 3 -27.8 5.59 9.0 -4.97 -40.4 -15.129 8e-04 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
>
> # using the esticon function from package doBy,
> # we can explicitly calculate our own lsmeans or whatever
> # this can be useful if we want/need to calculate our
> # own estimates of various things, and then use Scheffe's test
>
> w1 <- c(1,0,0)
> esticon(lm1,w1)
beta0 Estimate Std.Error X2.value DF Pr(>|X^2|) Lower Upper
1 0 333.4728 3.955878 7106.164 1 0 325.7195 341.2262
>
> w1 <- c(1,1,0)
> esticon(lm1,w1)
beta0 Estimate Std.Error X2.value DF Pr(>|X^2|) Lower Upper
1 0 346.3705 3.955878 7666.48 1 0 338.6171 354.1238
>
> w1 <- c(1,0,1)
> esticon(lm1,w1)
beta0 Estimate Std.Error X2.value DF Pr(>|X^2|) Lower Upper
1 0 374.1551 3.955878 8945.767 1 0 366.4017 381.9085
>
> # consider that we need to test the average of
> # treatments 1 and 2 vs trt 3. This is not provided
> # by either the lsmeans or the difflsmeans.
>
> # construct our own estimate/contrast
> w1 <- c(0,0.5,-1)
> diff1 <- esticon(lm1,w1)
> diff1
beta0 Estimate Std.Error X2.value DF Pr(>|X^2|) Lower Upper
1 0 -34.23344 4.844941 49.92569 1 0 -43.72935 -24.73753
>
> # get tabulated F value, df1 = ndf = Number of treatments - 1, = (3-1=2)
> # df2 = ddf = N trees - N treatments (12-3 = 9)
> Ftab <- qf(0.95,df1=2,df2=9)
> Ftab
[1] 4.256495
> ndf <- 2
> CD <- diff1$Std.Error * sqrt(ndf*Ftab)
> CD
[1] 14.1361
>
> # Estimate = -34.23, se = 4.8449
> # CD = 14.1361
> # |Estimate| > CD, therefore SS
> #
>
