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.padmanab...@itc.i
om> n>
cc
05/18/2 r-help@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.padmanab...@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.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
