[R] Query on linear mixed model

Joshua Wiley jwiley.psych at gmail.com
Tue May 18 15:48:46 CEST 2010


On Tue, May 18, 2010 at 4:52 AM, Vijayan Padmanabhan
<V.Padmanabhan at itc.in> wrote:
<snip>

This does not answer your statistical question, but I did include some
ideas to simplify your script.

> ##For ALL Product Comparison across All Time
> Points.
> options(contrasts=c('contr.treatment','contr.poly'))
> data<-subset(MyData,Attribute=="ChromaL")
> tapply(data$value, list(Product=data$Product),
> mean,
>  na.rm=TRUE)
> model <- lme(value ~
> Product+Time+Arm+Product*Arm+Product*Time+Product*Arm*Time,
>  random = ~1 | Subj,data =data)
> summary(model)

using '*' automatically crosses all the variables.  A more parsimonius form is:
lme(value ~ Product*Arm*Time, random = ~1 | Subj,data =data)

there is only a slight reordering of effects, but all estimates are the same.

> x<-anova(model)
> x
> library(multcomp)
> su<-summary(glht(model,linfct=mcp(Product="Tukey")))
> ##length(su)
> ##su[1:(length(su)-4)]
> x11()
> plot(summary(glht(model,linfct=mcp(Product="Tukey"))),cex.axis=0.6)
>
>
> ##For Each Product Comparison across All Time
> Points.
> data<-MyData
> data<-subset(data,Product=="a")
> tapply(data$value, list(Time=data$Time), mean,
>  na.rm=TRUE)
> model <- lme(value ~ Time+Arm+Time*Arm, random =
> ~1 | Subj,data =data)

again, simplified:

lme(value ~ Time*Arm, random = ~1 | Subj,data =data)

(no reordering even this time)

> summary(model)
> x<-anova(model)
> x
> library(multcomp)
> summary(glht(model,linfct=mcp(Time="Tukey")))
> x11()
> plot(summary(glht(model,linfct=mcp(Time="Tukey"))),cex.axis=0.6)
>
> data<-MyData
> data<-subset(data,Product=="b")
> tapply(data$value, list(Time=data$Time), mean,
>  na.rm=TRUE)
> model <- lme(value ~ Time+Arm+Time*Arm, random =
> ~1 | Subj,data =data)

lme(value ~ Time*Arm, random = ~1 | Subj,data =data)

> summary(model)
> x<-anova(model)
> x
> library(multcomp)
> summary(glht(model,linfct=mcp(Time="Tukey")))
> x11()
> plot(summary(glht(model,linfct=mcp(Time="Tukey"))),cex.axis=0.6)
>
> data<-MyData
> data<-subset(data,Product=="c")
> tapply(data$value, list(Time=data$Time), mean,
>  na.rm=TRUE)
> model <- lme(value ~ Time+Arm+Time*Arm, random =
> ~1 | Subj,data =data)

lme(value ~ Time*Arm, random = ~1 | Subj,data =data)

> summary(model)
> x<-anova(model)
> x
> library(multcomp)
> summary(glht(model,linfct=mcp(Time="Tukey")))
> x11()
> plot(summary(glht(model,linfct=mcp(Time="Tukey"))),cex.axis=0.6)
>
>
> ##For All Product Comparison at Each Time Points.
> data<-MyData
> data<-subset(data, Time=="T0")
> tapply(data$value, list(Product=data$Product),
> mean,
>  na.rm=TRUE)
>
> model <- lme(value ~ Product+Arm+Product:Arm,
> random = ~1 | Subj,data =data)

here you used ':' so it is not redundant, but it can still be simplified to:

lme(value ~ Product*Arm, random = ~1 | Subj,data =data)

> summary(model)
> x<-anova(model)
> x
> library(multcomp)
> summary(glht(model,linfct=mcp(Product="Tukey")))
> x11()
> plot(summary(glht(model,linfct=mcp(Product="Tukey"))),cex.axis=0.6)
>
>
> data<-MyData
> data<-subset(data, Time=="T1")
> tapply(data$value, list(Product=data$Product),
> mean,
>  na.rm=TRUE)
>
> model <- lme(value ~ Product+Arm+Product:Arm,
> random = ~1 | Subj,data =data)

lme(value ~ Product*Arm, random = ~1 | Subj,data =data)

> summary(model)
> x<-anova(model)
> x
> library(multcomp)
> summary(glht(model,linfct=mcp(Product="Tukey")))
> x11()
> plot(summary(glht(model,linfct=mcp(Product="Tukey"))),cex.axis=0.6)
>
>
> data<-MyData
> data<-subset(data, Time=="T2")
> tapply(data$value, list(Product=data$Product),
> mean,
>  na.rm=TRUE)
>
> model <- lme(value ~ Product+Arm+Product:Arm,
> random = ~1 | Subj,data =data)

lme(value ~ Product*Arm, random = ~1 | Subj,data =data)

> summary(model)
> x<-anova(model)
> x
> library(multcomp)
> summary(glht(model,linfct=mcp(Product="Tukey")))
> x11()
> plot(summary(glht(model,linfct=mcp(Product="Tukey"))),cex.axis=0.6)

Unless only parts of this script are being run in a given session, I
cannot think of a good reason to keep calling  "library(multcomp)".  I
also notice that your script repeats the same steps frequently, so you
might benefit from making a little function that does all the steps
you want.  That way if you ever want to add or change something, you
just have to update the function.

Good luck,

Josh

<snip>



-- 
Joshua Wiley
Senior in Psychology
University of California, Riverside
http://www.joshuawiley.com/



More information about the R-help mailing list