Nothing
#-*- 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())
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.