-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch]On Behalf Of Dominik
Grathwohl
Sent: Monday, March 07, 2005 17:18 PM
To: r-help at stat.math.ethz.ch
Subject: [R] One mixed effects model for two variables
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