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.padmanab...@itc.i r.edu> n> Sent cc by: r-help@r-project.org istazah Subject n...@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.padmanab...@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@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. ______________________________________________ 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.