[R] One mixed effects model for two variables

[Austin, Matt]

Here are two possibilities, but neither tests simultaneously--others may
have those suggestions.  

1.  If one of the variables in clinically more important than the other,
consider a stepdown procedure where the more important is tested first, then
if the first is successful the second is tested.  Procedures such as this
are well documented in the literature.

2.  Find a function of the two that has a clinical interpretation and model
that response. Body mass index is an example of such a function (kg/m^2),
but probably is not an endpoint that would be of interest in a growth study
in pediatrics.

Matt Austin
Statistician

Amgen 
One Amgen Center Drive
M/S 24-2-C
Thousand Oaks CA 93021
(805) 447 - 7431

"Today has the fatigue of a Friday and the desperation of a Monday"  -- S.
Pearce 2005


> -----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
> 
> ______________________________________________
> R-help at 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