[R] One mixed effects model for two variables

From: Dominik Grathwohl <DominikGrathwohl_at_gmx.de>
Date: Tue 08 Mar 2005 - 12:17:46 EST


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&#8217;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

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Tue Mar 08 12:32:57 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:30:41 EST