[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