69 lines
2.4 KiB
R
69 lines
2.4 KiB
R
#-*- R -*-
|
|
|
|
library( nlme )
|
|
options( width = 65, digits = 5 )
|
|
options( contrasts = c(unordered = "contr.helmert",
|
|
ordered = "contr.poly") )
|
|
pdf( file = 'ch02.pdf' )
|
|
|
|
# Chapter 2 Theory and Computational Methods for Linear Mixed-Effects Models
|
|
|
|
# 2.2 Likelihood Estimation for LME Models
|
|
|
|
Xmat <- matrix( c(1, 1, 1, 1, 8, 10, 12, 14), ncol = 2 )
|
|
Xmat
|
|
Xqr <- qr( Xmat ) # creates a QR structure
|
|
qr.R( Xqr ) # returns R
|
|
qr.Q( Xqr ) # returns Q-truncated
|
|
qr.Q( Xqr, complete = TRUE ) # returns the full Q
|
|
|
|
fm1Rail.lme <- lme( travel ~ 1, data = Rail, random = ~ 1 | Rail,
|
|
control = list( msVerbose = TRUE ) )
|
|
fm1Rail.lme <- lme( travel ~ 1, data = Rail, random = ~ 1 | Rail,
|
|
control = list( msVerbose = TRUE, niterEM = 0 ))
|
|
|
|
fm1Machine <-
|
|
lme( score ~ Machine, data = Machines, random = ~ 1 | Worker )
|
|
fm2Machine <- update( fm1Machine, random = ~ 1 | Worker/Machine )
|
|
anova( fm1Machine, fm2Machine )
|
|
|
|
OrthoFem <- Orthodont[ Orthodont$Sex == "Female", ]
|
|
fm1OrthF <- lme( distance ~ age, data = OrthoFem,
|
|
random = ~ 1 | Subject )
|
|
fm2OrthF <- update( fm1OrthF, random = ~ age | Subject )
|
|
orthLRTsim <- simulate.lme( fm1OrthF, m2 = fm2OrthF, nsim = 1000 )
|
|
plot( orthLRTsim, df = c(1, 2) ) # produces Figure 2.3
|
|
|
|
machineLRTsim <- simulate.lme(fm1Machine, m2 = fm2Machine, nsim= 1000)
|
|
plot( machineLRTsim, df = c(0, 1), # produces Figure 2.4
|
|
layout = c(4,1), between = list(x = c(0, 0.5, 0)) )
|
|
|
|
stoolLRTsim <-
|
|
simulate.lme( list(fixed = effort ~ 1, data = ergoStool,
|
|
random = ~ 1 | Subject),
|
|
m2 = list(fixed = effort ~ Type),
|
|
method = "ML", nsim = 1000 )
|
|
plot( stoolLRTsim, df = c(3, 4) ) # Figure 2.5
|
|
data( PBIB, package = 'SASmixed' )
|
|
pbibLRTsim <-
|
|
simulate.lme(list( fixed = response ~ 1, data = PBIB,
|
|
random = ~ 1 | Block ),
|
|
m2 = list(fixed = response ~ Treatment, data = PBIB,
|
|
random = ~ 1 | Block),
|
|
method = "ML", nsim = 1000 )
|
|
plot( pbibLRTsim, df = c(14,16,18), weights = FALSE ) # Figure 2.6
|
|
|
|
summary( fm2Machine )
|
|
|
|
fm1PBIB <- lme(response ~ Treatment, data = PBIB, random = ~ 1 | Block)
|
|
anova( fm1PBIB )
|
|
fm2PBIB <- update( fm1PBIB, method = "ML" )
|
|
fm3PBIB <- update( fm2PBIB, response ~ 1 )
|
|
anova( fm2PBIB, fm3PBIB )
|
|
anova( fm2Machine )
|
|
|
|
# cleanup
|
|
|
|
summary(warnings())
|
|
|