[R] Query on linear mixed model

[Ista Zahn]

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