[R] Query on linear mixed model

Vijayan Padmanabhan V.Padmanabhan at itc.in
Wed May 19 05:51:23 CEST 2010


Thanks Joshua..
It really helped in polishing my coding essentials
in R.
Thanks & Regards
Vijayan Padmanabhan


"What is expressed without proof can be denied
without proof" - Euclide.



                                                  
             Joshua                               
             Wiley                                
             <jwiley                           To 
             .psych@         Vijayan Padmanabhan  
             gmail.c         <V.Padmanabhan at itc.i 
             om>             n>                   
                                               cc 
             05/18/2         r-help at r-project.org 
             010                          Subject 
             07:18           Re: [R] Query on     
             PM              linear mixed model   
                                                  
                                                  
                                                  
                                                  
                                                  
                                                  




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/



Can you avoid printing this?
Think of the environment before printing the email.
-------------------------------------------------------------------------------
Please visit us at www.itcportal.com
******************************************************************************
This Communication is for the exclusive use of the intended recipient (s) and shall
not attach any liability on the originator or ITC Ltd./its Subsidiaries/its Group 
Companies. If you are the addressee, the contents of this email are intended for your 
use only and it shall not be forwarded to any third party, without first obtaining 
written authorisation from the originator or ITC Ltd./its Subsidiaries/its Group 
Companies. It may contain information which is confidential and legally privileged
and the same shall not be used or dealt with by any third party in any manner 
whatsoever without the specific consent of ITC Ltd./its Subsidiaries/its Group 
Companies.



More information about the R-help mailing list