[R] Query on linear mixed model

[Vijayan Padmanabhan]

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.