example: lmm package example command file

Description See Also Examples

Description

The data as contained in marijuana is used to fit a compound symmetry model with a fixed effect for each occasion and a random intercept for each subject.

Since the six measurements per subject were not clearly ordered in time, instead of a model with time of measurement entered with linear (or perhaps higher-order polynomial) effects, the model has an intercept and five dummy codes to allow the population means for the six occasions to be estimated freely. For a subject i with no missing values, the covariate matrices will be

Xi=( 1 1 0 0 0 0, 1 0 1 0 0 0, 1 0 0 1 0 0, 1 0 0 0 1 0, 1 0 0 0 0 1) and Zi=(1,1,1,1,1,1)

The Xi's and Zi's are combined into a single matrix called pred (Zi is merely the first column of Xi), simply the matrices Xi (i=1,...,9), stacked upon each other.

See Also

ecmeml, ecmerml, fastml, fastrml, fastmcmc, fastmode, mgibbs

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
### Model specification ###
data(marijuana)
# To work only on those with complete data
marijuana <- subset(marijuana,!is.na(y))
attach(marijuana)
pred <- cbind(int,dummy1,dummy2,dummy3,dummy4,dummy5)
xcol <- 1:6
zcol <- 1

### ML Estimation ###
ecmeml.result <- ecmeml(y,subj,pred,xcol,zcol)
fastml.result <- fastml(y,subj,pred,xcol,zcol)
#
# which converged in 212 and 8 cycles, respectively. For example, the
# first elemenent of the ML estimate of the fixed effects (the intercept)
# estimates the mean for the last occasion and the other elements of beta
# estimate the differences in means between the first five occasions and
# the last one. So we can find the estimated means for the six occasions.
#
beta.hat <- fastml.result$beta
muhat <- c(beta.hat[2]+beta.hat[1], beta.hat[3]+beta.hat[1],
   beta.hat[4]+beta.hat[1], beta.hat[5]+beta.hat[1],
   beta.hat[6]+beta.hat[1], beta.hat[1]) 

### RML estimation ###
ecmerml.result <- ecmerml(y,subj,pred,xcol,zcol)
fastrml.result <- fastrml(y,subj,pred,xcol,zcol)

### Improved variance estimation in Section 4 ### 
b.hat <- as.vector(fastrml.result$b.hat)
se.new <- sqrt(as.vector(fastrml.result$cov.b.new))
se.old <- sqrt(as.vector(fastrml.result$cov.b))
table2 <- cbind(round(b.hat,3),round(cbind(b.hat-2*se.old,b.hat+2*se.old,
      b.hat-2*se.new,b.hat+2*se.new),2),round(100*(se.new-se.old)/se.old))
dimnames(table2) <- list(paste("Subject",format(1:9)),
   c("Est.","Lower.old","Upper.old","Lower.new","Upper.new","Increase (%)"))
print(table2)
#
# which reproduces Table 2 and compares 95% interval estimates
# under the new method to conventional empirical Bayes intervals.

### MCMC in Section 5 ###
prior <- list(a=3*100,b=3,c=3,Dinv=3*5)
gibbs.result <- mgibbs(y,subj,pred,xcol,zcol,prior=prior,seed=1234,iter=5000)
fmcmc.result <- fastmcmc(y,subj,pred,xcol,zcol,prior=prior,seed=2345,iter=5000)
#
# which run 5,000 cycles for each algorithm and generates Figure 1.
#
# library(ts)
par(mfrow=c(2,1))
acf(log(gibbs.result$psi.series[1,1,]),lag.max=10, ylim=0:1)
acf(log(fmcmc.result$psi.series[1,1,]),lag.max=10, ylim=0:1)
detach(marijuana)

Example output

Performing ECME...
Performing FAST-ML...
Performing ECME for RML estimation...
Performing FAST-RML...
            Est. Lower.old Upper.old Lower.new Upper.new Increase (%)
Subject 1  0.190     -2.39      2.77     -3.21      3.59           32
Subject 2 -0.168     -2.75      2.41     -3.43      3.09           26
Subject 3  0.011     -2.57      2.59     -2.59      2.61            1
Subject 4  0.215     -2.40      2.83     -3.46      3.89           40
Subject 5 -0.459     -3.06      2.14     -6.50      5.58          133
Subject 6 -0.288     -2.87      2.29     -4.54      3.96           65
Subject 7  0.668     -1.91      3.25     -7.52      8.86          218
Subject 8 -0.482     -3.06      2.10     -6.68      5.71          140
Subject 9  0.313     -2.31      2.93     -4.27      4.90           75
Performing modified Gibbs...
Performing FAST-MODE...
Performing FAST-MCMC...

lmm documentation built on July 8, 2020, 6:28 p.m.