library(nlme)
ex.data<-read.table(file="W:/teaching/stat579/data/oncology2.dat",head=T)
#ex.data<-read.table(file="C:/jenn/teachingrelated/stat579/mixedmodel1/oncology2.dat",head=T)
sort.data<-ex.data[order(ex.data$oncologist, ex.data$shape), ]
#define global variables
oncologist<-as.factor(sort.data$oncologist)
shape<-as.factor(sort.data$shape)
material<-as.factor(sort.data$material)
response<-sort.data$response
#ex.data<-read.table(file="W:/teaching/stat579/data/oncology2.dat",head=T,colClasses=c('numeric','factor','factor',
#'factor','numeric'))
nrow(sort.data)
head(sort.data)
#calculate mean value of response of three shapes
data1<-sort.data[which(sort.data$shape==1),]
#colClasses=c('numeric','factor','factor','factor','factor')]
mean1<-mean(data1$response)
mean1
data2<-sort.data[which(sort.data$shape==2),]
mean2<-mean(data2$response)
mean2
data3<-sort.data[which(sort.data$shape==3),]
mean3<-mean(data3$response)
mean3

#graphs
par(mfrow=c(2,2))
boxplot(response~shape,main="shape vs. Response", xlab="shape",ylab="response")
boxplot(response~material,main="material vs. Response", xlab="material",ylab="response")
boxplot(response~oncologist,main="oncologist vs. Response", xlab="oncologist",ylab="response")
library(lattice)
trellis.device(color=F)
xyplot(response ~ oncologist, group=shape, main="profile plot")

#start with full model, check if interaction term can be deleted
library(lme4)
model1=lmer(response~shape+(1|oncologist)+(1|shape:oncologist)) 
model2=lmer(response~shape+(1|oncologist)) 
anova(model1,model2)

#reduced model
inf.lme<-lme(response~shape,random=~1|oncologist)
summary(inf.lme)
fixef(inf.lme)
ranef(inf.lme)
fit.lme<-predict(inf.lme)
resid.lme<-resid(inf.lme)
par(mfrow=c(2,2))
plot(resid.lme~fit.lme,main="fitted values vs. residuals",xlab="fitted",ylab="residuals")
abline(0,0,lty=2)
qqnorm(resid.lme)
qqline(resid.lme)
hist(resid.lme,main="histogram of residuals",xlab="residuals")

inf.lme2<-gls(response~shape)
anova(inf.lme,inf.lme2)