[R] One mixed effects model for two variables

[Dominik Grathwohl]

Hello R-help group,

I need some suggestions in stating a mixed effects model. I would like to
fit one mixed effect model to two or more variables and than investigating
treatment effects by applying the multcomp library. I will explain my
problem along an example of weight and length measurements on infants.
Weight and length  are measured at two occasions in time and for two
treatment groups, one group got some experimental formula and the other some
control formula. It is also reasonable to assume that weight and length are
correlated. Aim is to estimate the treatment effect on weight and on length
taking multiplicity into account. 
Below I generated some data of the proposed properties and I applied also
two mixed effects models, one model for weight and one for length. Up to now
I’m not able to state a model for both variables together, any
suggestions?

library(nlme)

n <- 200                # n = number of subjects
id <- rep(1:n, each=2)  # id = subject identification number
t <- rep(0:1, times=n)  # two occations in time
trt <- rep(sample(rep(0:1, times=n/2)), each=2) # two treatment groups
b1 <- rnorm(n, mean=0, sd=0.5)              # b1 = random effect of weight

y1 <- 3.5 + rep(b1, each=2) + 0.7 * t + 0.01 * trt + 0.1 * t * trt +
rnorm(2*n,mean=0, sd=0.09)
# y1 = weight measurements, 3.5 kg at baseline, 
# 0.7 kg more at the second time occation 
# 0.01 = the treatment effect at baseline, 
# the treatment effect is 0.1 kg for the 
# experimental group at at the second time occation. 
# The whithin subject standard deviation is 0.09.

b2 <- 0.9 * 3/sd(b1) * b1 + rnorm(length(b1), mean=0, sd=sqrt(1-0.9**2)*3)
# b2 = random effect for length with standard deviation = 3 
# and correlated with the random effect of weight (b1) by r=0.9

y2 <- 49 + rep(b2, each=2) + 2 * t - 0.05 * trt + 0.5 * t * trt +
rnorm(2*n,mean=0, sd=1)
# y2 = length measurements, 49 cm at baseline, 
# 2 cm more at the second time occation 
# 0.05 = the treatment effect at baseline, 
# the treatment effect is 0.5 cm for the 
# experimental group at at the second time occation. 
# The whithin subject standard deviation is 1.


# data frame of the data:
df <- data.frame(var=as.factor(rep(0:1, each=2*n)),
 id=c(id, id), t=as.factor(c(t, t)), trt=as.factor(c(trt, trt)), y=c(y1,
y2))

# grouped data object:
gd <- groupedData(y ~ t | id, data=df)

# mixed effects model on weight:
fm1weight <- lme(y ~ t * trt, random = ~ 1 | id, data = gd, subset=var==0)
summary(fm1weight)

# mixed effect model on length:
fm1length <- lme(y ~ t * trt, random = ~ 1 | id, data = gd, subset=var==1)
summary(fm1length)

-- 
DSL Komplett von GMX +++ Supergünstig und stressfrei einsteigen!
AKTION "Kein Einrichtungspreis" nutzen: http://www.gmx.net/de/go/dsl