[R] Query on linear mixed model

Wed May 19 05:52:33 CEST 2010

Thanks Ista..
I will take your suggestion.
Regards
Vijayan Padmanabhan

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

             Ista                                 
             Zahn                                 
             <izahn@                           To 
             psych.r         Vijayan Padmanabhan  
             ocheste         <V.Padmanabhan at itc.i 
             r.edu>          n>                   
             Sent                              cc 
             by:             r-help at r-project.org 
             istazah                      Subject 
             n at gmail         Re: [R] Query on     
             .com            linear mixed model   

             05/18/2                              
             010                                  
             06:39                                
             PM                                   

Hi Vijayan,
You are really asking for this list to serve as
your statistical
consultant, which is not its purpose. If you have
a specific problem
(and if you know how to ask for help -- see the
posting guide) this
list is a tremendous resource. But it is not a
replacement for a
statistician.

Best,
Ista

On Tue, May 18, 2010 at 7:52 AM, Vijayan
Padmanabhan
<V.Padmanabhan at itc.in> wrote:
>
>
> Hi R Forum
> I am a newbie to R and I have been amazed by
what
> I can get my team to accomplish just by
> implementing Scripting routines of R in all my
> team's areas of interest..
> Recently i have been trying to adopt R scripting
> routine for some analysis with longitudanal
data..
> I am presenting my R script below that I have
> tried to make to automate data analysis for
> longitudanal data by employing functions from
> library(nlme) and library(multcomp)..
>
> I would be thankful for receiving inputs on this
> script and let me know if I have modeled the lme
> formula correctly.. If the formula i have used
is
> not the correct one i would appreciate receiving
> inputs on what is the correct formula for lme
that
> I should use given the context of this example
> study shown in the data..
>
> Just to give an introduction.. the data is about
> studying efficacy of 3 products for their impact
> in skin brightness over 3 time points of
> investigation .. The design is such that all the
> products are tried on patches made on each arm
> (left and right) for each volunteer and chromaL
is
> measured as the response over 3 time points
> Baseline (referred as T0), T1 and T2..
>
>
> The answers i am looking to get from the
analysis
> routine is as follows:
>
> Overall across different time points studied
which
> products is superior?
> For Each Product is their a significant
difference
> in the response variable across different time
> points of investigation?
> For Each Time Point is their a significant
> difference between the different products for
the
> measured response?
> Regards
> Vijayan Padmanabhan
> Research Scientist, ITC R&D, Phase I, Peenya,
> Bangalore - 58
>
> The Full R script is given below:
>
> MyData <- data.frame(Subj = factor(rep(c("S1",
> "S2", "S3"), 18)),
> Product = factor(rep(letters[1:3],each=3,54)),
> Arm = factor(rep(c("L","R"),each=9,54)),
> Attribute =
factor(rep(c("ChromaL"),each=54,54)),
> Time =
factor(rep(c("T0","T1","T2"),each=18,54)),
> value=as.numeric(c(43.52,
>
44.22,43.2,40.12,39.45,38.98,43.23,42.34,44.54,50.23,45.12,43.56,43.23,44.56,42.34,45.67,

>
43.23,44.54,43.52,44.22,43.2,40.12,39.45,38.98,43.23,42.34,44.54,50.23,45.12,43.56,

>  43.23, 44.56, 42.34, 45.67, 43.23,
> 44.54, 45.5, 46.45, 47.56, 46.98, 46.3, 43.1,
> 45.6, 44.2, 40.1, 49.8, 48, 46, 47.98, 46.9,
> 43.78, 47.35, 44.9, 48)))
> tapply(MyData$value,
> list(Attribute=MyData$Attribute,
>  Product=MyData$Product), mean, na.rm=TRUE)
>
>
> Time = factor(MyData$Time)
> Product = factor(MyData$Product)
> Subj = factor(MyData$Subj)
> Attribute=factor(MyData$Attribute)
> Arm=factor(MyData$Arm)
> ##library(reshape)
> ##data<-melt(data, id=c("Product",
> "Subject","Attribute"))
> ##data$Product<-as.character(Data$Product)
> library(nlme)
> library(multcomp)
>
>
> ##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)
> 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)
> 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)
> 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)
> 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)
> 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)
> 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)
> 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)

>
>
>
###################################################################################

> ##Regards
> ##Vijayan Padmanabhan
> Research Scientist, ITC R&D, Bangalore-58
> "What is expressed without proof can be denied
> without proof" - Euclide.
>
>
> 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 at 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.
>

--
Ista Zahn
Graduate student
University of Rochester
Department of Clinical and Social Psychology
http://yourpsyche.org

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.