#############################################
# R Codes for paper "Modelling finger force produced from different tasks 
# using linear mixed models with lme R function"
# Caroline Bazzoli, Frédérique Letué, Marie-José Martinez
# Université de Grenoble, LJK
# Conforme Latex, 27/05/2014
#############################################
library(nlme)
###################
# 2. The data     #
###################

#loading the data
read.table("data.txt",header=TRUE)-> Data
attach(Data)
summary(Data)

# Defining the variables
F<-c(I.ExtP3,M.ExtP3,R.ExtP3,L.ExtP3,I.FlexP3,M.FlexP3,R.FlexP3,
     L.FlexP3,I.ExtP1,M.ExtP1,R.ExtP1,L.ExtP1)
location <- gl(3,180,540,labels=c("ExtP3","FlexP3","ExtP1"))
finger <- gl(4,45,540,labels=c("I","M","R","L"))
individual <- as.factor(rep(rep(1:15,each=3),12))
trial <- as.factor(c(rep(1:45,4),rep(1:45,4)+45,rep(1:45,4)+90))
Data.new<-data.frame(F,location,finger,individual,trial)
head(Data.new,200)
###############################
# 3.Preliminary study         #
#################################
# 3.1 Exploratory data analysis #
#################################
# raw data by location, finger and individual
# ***Figure 1***
op<-par(mfrow=c(1,3))
plot(I.ExtP3~as.numeric(individual)[(location=="ExtP3")&(finger=="I")],
     col="blue",pch=1,
     ylim=c(0,50),main="Location ExtP3",
     xlab="individual",ylab="Finger force intensity")
points(M.ExtP3~individual[(location=="ExtP3")&(finger=="M")],
       col="red",pch=2)
points(R.ExtP3~individual[(location=="ExtP3")&(finger=="R")],
       col="green",pch=3)
points(L.ExtP3~individual[(location=="ExtP3")&(finger=="L")],
       col="magenta",pch=4)
legend("topleft",c("I","M","R","L"),pch=1:4,
       col=c("blue","red","green","magenta"))

plot(I.FlexP3~as.numeric(individual)[(location=="FlexP3")&(finger=="I")],
     col="blue",pch=1,
     ylim=c(0,50),main="Location FlexP3",
     xlab="individual",ylab="Finger force intensity")
points(M.FlexP3~individual[(location=="ExtP3")&(finger=="M")],
       col="red",pch=2)
points(R.FlexP3~individual[(location=="ExtP3")&(finger=="R")],
       col="green",pch=3)
points(L.FlexP3~individual[(location=="ExtP3")&(finger=="L")],
       col="magenta",pch=4)
legend("topleft",c("I","M","R","L"),pch=1:4,
       col=c("blue","red","green","magenta"))

plot(I.ExtP1~as.numeric(individual)[(location=="ExtP1")&(finger=="I")],
     col="blue",pch=1,
     ylim=c(0,50),main="Location ExtP1",
     xlab="individual",ylab="Finger force intensity")
points(M.ExtP1~individual[(location=="ExtP3")&(finger=="M")],
       col="red",pch=2)
points(R.ExtP1~individual[(location=="ExtP3")&(finger=="R")],
       col="green",pch=3)
points(L.ExtP1~individual[(location=="ExtP3")&(finger=="L")],
       col="magenta",pch=4)
legend("topleft",c("I","M","R","L"),pch=1:4,
       col=c("blue","red","green","magenta"))
#***
# finger correlation 
#***Figure 2***
op<-par(mfrow=c(1,1))
index<-c(I.ExtP3,I.FlexP3,I.ExtP1)
middle<-c(M.ExtP3,M.FlexP3,M.ExtP1)
ring<-c(R.ExtP3,R.FlexP3,R.ExtP1)
little<-c(L.ExtP3,L.FlexP3,L.ExtP1)
cor(cbind(index,middle,ring,little))
plot(data.frame(index,middle,ring,little),pch=rep(1:3,each=60),
     col=rep(1:3,each=60))
#***
#######################################
# 3.2 Two-factor ANOVA and its limits #
#######################################
fit<-lm(F~location*finger)
anova(fit)
summary(fit)
plot(fit)

# residuals
#***Figure 3***
res <- fit$residuals
col <- as.numeric(finger)
pch <- as.numeric(finger)
op<-par(mfrow=c(1,3))
plot(res[location=="ExtP3"]~as.numeric(individual)[location=="ExtP3"],
     col=col,pch=pch,
     ylim=c(-21,21), main="Location ExtP3",
     xlab="Individual",ylab="residuals")
legend("topright",c("I","M","R","L"),pch=1:4,
       col=c("blue","red","green","magenta"))

plot(res[location=="FlexP3"]~as.numeric(individual)[location=="FlexP3"],
     col=col,pch=pch,
     ylim=c(-21,21), main="Location FlexP3",
     xlab="Individual",ylab="residuals")
legend("topleft",c("I","M","R","L"),pch=1:4,
       col=c("blue","red","green","magenta"))

plot(res[location=="ExtP1"]~as.numeric(individual)[location=="ExtP1"],
     col=col,pch=pch,
     ylim=c(-21,21),main="Location ExtP1",
     xlab="Individual",ylab="residuals")
legend("topright",c("I","M","R","L"),pch=1:4,
       col=c("blue","red","green","magenta"))
#***
# residuals correlation
#***Figure 4***
op<-par(mfrow=c(1,1))
index.res<-res[finger=="I"]
middle.res<-res[finger=="M"]
ring.res<-res[finger=="R"]
little.res<-res[finger=="L"]
cor(cbind(index.res,middle.res,ring.res,little.res))
plot(data.frame(index.res,middle.res,ring.res,little.res),
     pch=rep(1:3,each=60),col=rep(1:3,each=60))
#***
#############################################################
# 4. Model specification using a linear mixed-effects model #
#############################################################
# 4.1 Modelling the random effect structure #
#############################################
#***Table 1***
library(nlme)
fitM0 <- lme(F ~ finger*location, random=~1|individual, method="ML")
summary(fitM0)
resM0.std<-residuals(fitM0,type="pearson")
#***Figure 5***
plot(fitM0,individual~resM0.std|location,abline=0,xlim=c(-5,5),
     xlab="Standardized residuals")
#***Figure 6***
plot(fitM0,individual~resM0.std|finger,abline=0,xlim=c(-5,5),
     xlab="Standardized residuals")
#***
# individual random effects with identical variances
#***Table 2***
fitM1 <- lme(F ~ finger*location, 
            random=list(individual=pdBlocked(list(pdIdent(~1),
                                                  pdIdent(~location-1),
                                                  pdIdent(~finger-1),
                                                  pdIdent(~location:finger-1)))),
            method="ML")
resM1.std<-residuals(fitM1,type="pearson")
#***Figure 7***
plot(fitM1,individual~resM1.std|location,abline=0,xlim=c(-5,5),
     xlab="Standardized residuals")
#***Figure 8***
plot(fitM1,individual~resM1.std|finger,abline=0,xlim=c(-5,5),
     xlab="Standardized residuals")
#***
#***Table 3***
fitM2 <- lme(F ~ finger*location, 
               random=list(individual=pdBlocked(list(pdIdent(~1),
                                                     pdDiag(~location-1),
                                                     pdIdent(~finger-1),
                                                     pdIdent(~location:finger-1)))),
               method="ML",control=lmeControl(msMaxIter=1000))

resM2.std<-residuals(fitM2,type="pearson")
#***Table 4***
anova(fitM0,fitM1,fitM2)
# residuals correlation from model M2
#***Figure 9***
op<-par(mfrow=c(1,1))
index.resM2.std<-resM2.std[finger=="I"]
middle.resM2.std<-resM2.std[finger=="M"]
ring.resM2.std<-resM2.std[finger=="R"]
little.resM2.std<-resM2.std[finger=="L"]
cor(cbind(index.resM2.std,middle.resM2.std,ring.resM2.std,little.resM2.std))
plot(data.frame(index.resM2.std,middle.resM2.std,ring.resM2.std,little.resM2.std),
     pch=rep(1:3,each=60),col=rep(1:3,each=60))
#***
############################################################
# 4.2 Modelling the residual variance-covariance structure #
############################################################
# 4.2.1  Modelling the variance matrix for each location #
##########################################################
#***Table 5***
fitM2.1 <- lme(F ~ finger*location, 
             random=list(individual=pdBlocked(list(pdIdent(~1),
                                                   pdDiag(~location-1),
                                                   pdIdent(~finger-1),
                                                   pdIdent(~location:finger-1)))),
             weights=varIdent(form=~1|location),
             method="ML",control=lmeControl(msMaxIter=1000))
#***Table 6***
anova(fitM2,fitM2.1)
resM2.1.std<-residuals(fitM2.1,type="pearson")
#***Figure 10***
op<-par(mfrow=c(1,2))
boxplot(resM2.std~location,xlab="Location",
        ylab="Standardized residuals",
        main="Standardized residuals vs Location for M2",
        ylim=c(-4,4))
boxplot(resM2.std~finger,xlab="Finger",
        ylab="Standardized residuals",
        main="Standardized residuals vs Finger for M2",
        ylim=c(-4,4))
#***Figure 11***
op<-par(mfrow=c(1,2))
boxplot(resM2.1.std~location,xlab="Location",
        ylab="Standardized residuals",
        main="Standardized residuals vs Location for M2.1",
        ylim=c(-4,4))
boxplot(resM2.1.std~finger,xlab="Finger",
        ylab="Standardized residuals",
        main="Standardized residuals vs Finger for M2.1",
        ylim=c(-4,4))
#***
Index<-(finger=="I")
Index[Index==TRUE]<-"I"
Index[Index==FALSE]<-"other"
fitM2.2 <- lme(F ~ finger*location, 
               random=list(individual=pdBlocked(list(pdIdent(~1),
                                                     pdDiag(~location-1),
                                                     pdIdent(~finger-1),
                                                     pdIdent(~location:finger-1)))),
               weights=varIdent(form=~1|location*Index),
               method="ML",control=lmeControl(msMaxIter=1000))
anova(fitM2.1,fitM2.2)
resM2.2.std<-residuals(fitM2.2,type="pearson")
#***Figure 12***
op<-par(mfrow=c(1,2))
boxplot(resM2.2.std~location,xlab="Location",
        ylab="Standardized residuals",
        main="Standardized residuals vs Location for M2.2",
        ylim=c(-4,4))
boxplot(resM2.2.std~finger,xlab="Finger",
        ylab="Standardized residuals",
        main="Standardized residuals vs Finger for M2.2",
        ylim=c(-4,4))
#***
# residuals correlation from model M2.2
#***Table 7***
op<-par(mfrow=c(1,1))
index.resM2.2.std<-resM2.2.std[finger=="I"]
middle.resM2.2.std<-resM2.2.std[finger=="M"]
ring.resM2.2.std<-resM2.2.std[finger=="R"]
little.resM2.2.std<-resM2.2.std[finger=="L"]
cor(cbind(index.resM2.2.std,middle.resM2.2.std,ring.resM2.2.std,little.resM2.2.std))
plot(data.frame(index.resM2.2.std,middle.resM2.2.std,ring.resM2.2.std,little.resM2.2.std),
     pch=rep(1:3,each=60),
     main="Residuals correlation by finger",col=rep(1:3,each=60))
#***

#############################################
# 4.2.2  Modelling the correlation matrix   #
#############################################
fitM2.3 <- lme(F ~ finger*location, 
                         random=list(individual=pdBlocked(list(pdIdent(~1),
                                                               pdDiag(~location-1),
                                                               pdIdent(~finger-1),
                                                               pdIdent(~location:finger-1)))),                             
                         weights=varIdent(form=~1|location*Index),
                         correlation=corSymm(form=~1|individual/trial),
                         method="ML", control=lmeControl(msMaxIter=1000))
#***Table 8***
anova(fitM2.2, fitM2.3)

resM2.3.norm<-residuals(fitM2.3,type="normalized")
#***Figure 13***
op<-par(mfrow=c(1,2))
boxplot(resM2.3.norm~location,xlab="Location",
        ylab="Normalized residuals",
        main="Normalized residuals vs Location for M2.3")
boxplot(resM2.3.norm~finger,xlab="Finger",
        ylab="Normalized residuals",
        main="Normalized residuals vs Finger for M2.3")
#***
# residuals correlation from model M2.3
#***Table 9***
op<-par(mfrow=c(1,1))
index.resM2.3.norm<-resM2.3.norm[finger=="I"]
middle.resM2.3.norm<-resM2.3.norm[finger=="M"]
ring.resM2.3.norm<-resM2.3.norm[finger=="R"]
little.resM2.3.norm<-resM2.3.norm[finger=="L"]
cor(cbind(index.resM2.3.norm,middle.resM2.3.norm,ring.resM2.3.norm,little.resM2.3.norm))
plot(data.frame(index.resM2.3.norm,middle.resM2.3.norm,ring.resM2.3.norm,little.resM2.3.norm),
     pch=rep(1:3,each=60),
     main="Residuals correlation by finger from model M2.3",col=rep(1:3,each=60))
#***
# finger correlation location by location  
#***Table 10***
index.ExtP3<-resM2.3.norm[finger=="I"&location=="ExtP3"]
middle.ExtP3<-resM2.3.norm[finger=="M"&location=="ExtP3"]
ring.ExtP3<-resM2.3.norm[finger=="R"&location=="ExtP3"]
little.ExtP3<-resM2.3.norm[finger=="L"&location=="ExtP3"]

index.FlexP3<-resM2.3.norm[finger=="I"&location=="FlexP3"]
middle.FlexP3<-resM2.3.norm[finger=="M"&location=="FlexP3"]
ring.FlexP3<-resM2.3.norm[finger=="R"&location=="FlexP3"]
little.FlexP3<-resM2.3.norm[finger=="L"&location=="FlexP3"]

index.ExtP1<-resM2.3.norm[finger=="I"&location=="ExtP1"]
middle.ExtP1<-resM2.3.norm[finger=="M"&location=="ExtP1"]
ring.ExtP1<-resM2.3.norm[finger=="R"&location=="ExtP1"]
little.ExtP1<-resM2.3.norm[finger=="L"&location=="ExtP1"]

print("ExtP3")
cor(cbind(index.ExtP3,middle.ExtP3,ring.ExtP3,little.ExtP3))
print("FlexP3")
cor(cbind(index.FlexP3,middle.FlexP3,ring.FlexP3,little.FlexP3))
print("ExtP1")
cor(cbind(index.ExtP1,middle.ExtP1,ring.ExtP1,little.ExtP1))
#***
############################################################
# 5  Results  #
############################################################
fitM2.3.REML <- lme(F ~ finger*location, 
                    random=list(individual=pdBlocked(list(pdIdent(~1),
                                                          pdDiag(~location-1),
                                                          pdIdent(~finger-1),
                                                          pdIdent(~location:finger-1)))),                             
                    weights=varIdent(form=~1|location*Index),
                    correlation=corSymm(form=~1|individual/trial),
                    method="REML", control=lmeControl(msMaxIter=1000))
#############################################
# 5.1  Residual analysis of the final model #
#############################################
resM2.3.REML<-residuals(fitM2.3.REML,type="normalized")
#***Figure 14***
op<-par(mfrow=c(2,2))
hist(resM2.3.REML,freq=FALSE,xlab="Normalized residuals",
     main="Histogram of the normalized residuals")
lines((-400:400)/100,dnorm((-400:400)/100))
qqnorm(resM2.3.REML)            
qqline(resM2.3.REML)            
plot(resM2.3.REML~fitted(fitM2.3.REML),xlab="Fitted values",
     ylab="Normalized residuals",
     main="Normalized residuals vs fitted values")
plot(resM2.3.REML~F,xlab="Observed values",
     ylab="Normalized residuals",
     main="Normalized residuals vs observed values")
#***
#########################
# 5.2  Results analysis #
#########################
#***Table 11***
summary(fitM2.3.REML)
# Residual standard deviation
#***Table 12***
fitM2.3.REML$sigma*unique(1/varWeights(fitM2.3.REML$modelStruct$varStruct))
# estimated means levels of the finger-location crossing groups
#***Table 13***
fitM2.3.REML$coefficients$fixed
fitM2.3.REML$coefficients$fixed[1] # ExtP3:I
fitM2.3.REML$coefficients$fixed[1]+fitM2.3.REML$coefficients$fixed[2] #ExtP3:M
fitM2.3.REML$coefficients$fixed[1]+fitM2.3.REML$coefficients$fixed[3] #ExtP3:R
fitM2.3.REML$coefficients$fixed[1]+fitM2.3.REML$coefficients$fixed[4] #ExtP3:L

fitM2.3.REML$coefficients$fixed[1]+fitM2.3.REML$coefficients$fixed[5] # FlexP3:I
fitM2.3.REML$coefficients$fixed[1]+fitM2.3.REML$coefficients$fixed[2]+
  fitM2.3.REML$coefficients$fixed[5] +fitM2.3.REML$coefficients$fixed[7]#FlexP3:M
fitM2.3.REML$coefficients$fixed[1]+fitM2.3.REML$coefficients$fixed[3]+
  fitM2.3.REML$coefficients$fixed[5] +fitM2.3.REML$coefficients$fixed[8]#FlexP3:R
fitM2.3.REML$coefficients$fixed[1]+fitM2.3.REML$coefficients$fixed[4]+
  fitM2.3.REML$coefficients$fixed[5] + fitM2.3.REML$coefficients$fixed[9]#FlexP3:L

fitM2.3.REML$coefficients$fixed[1]+fitM2.3.REML$coefficients$fixed[6] # ExtP1:I
fitM2.3.REML$coefficients$fixed[1]+fitM2.3.REML$coefficients$fixed[2]+
  fitM2.3.REML$coefficients$fixed[6] +fitM2.3.REML$coefficients$fixed[10]# ExtP1:M
fitM2.3.REML$coefficients$fixed[1]+fitM2.3.REML$coefficients$fixed[3]+
  fitM2.3.REML$coefficients$fixed[6] +fitM2.3.REML$coefficients$fixed[11]# ExtP1:R
fitM2.3.REML$coefficients$fixed[1]+fitM2.3.REML$coefficients$fixed[4]+
  fitM2.3.REML$coefficients$fixed[6] + fitM2.3.REML$coefficients$fixed[12] # ExtP1:L

# contrastes
#***Table 14***
group <- gl(12,45,540,labels=c("ExtP3:I","ExtP3:M","ExtP3:R","ExtP3:L",
                               "FlexP3:I","FlexP3:M","FlexP3:R","FlexP3:L",
                               "ExtP1:I","ExtP1:M","ExtP1:R","ExtP1:L"))

library(MASS)
M.location<-cbind(
  c(1,0,0,0,-1,0,0,0,0,0,0,0), # ExtP3/FlexP3,I
  c(0,1,0,0,0,-1,0,0,0,0,0,0), # ExtP3/FlexP3,M
  c(0,0,1,0,0,0,-1,0,0,0,0,0), # ExtP3/FlexP3,R
  c(0,0,0,1,0,0,0,-1,0,0,0,0), # ExtP3/FlexP3,L
  c(1,0,0,0,0,0,0,0,-1,0,0,0), # ExtP3/ExtP1,I
  c(0,1,0,0,0,0,0,0,0,-1,0,0), # ExtP3/ExtP1,M
  c(0,0,1,0,0,0,0,0,0,0,-1,0), # ExtP3/ExtP1,R
  c(0,0,0,1,0,0,0,0,0,0,0,-1) # ExtP3/ExtP1,L
  )
contrasts(group)<-t(ginv(M.location))
fitM2.3.REML.location <- lme(F ~ group,
random=list(individual=pdBlocked(list(pdIdent(~1),
pdDiag(~location-1),
pdIdent(~finger-1),
pdIdent(~location:finger-1)))),
weights=varIdent(form=~1|location*Index),
correlation=corSymm(form=~1|individual/trial),
method="REML", control=lmeControl(msMaxIter=1000))
#***Table 15***
summary(fitM2.3.REML.location)

M.finger<-cbind(
  c(1,-1,0,0,0,0,0,0,0,0,0,0), # I/M, ExtP3
  c(0,0,0,0,1,-1,0,0,0,0,0,0), # I/M, FlexP3
  c(0,0,0,0,0,0,0,0,1,-1,0,0), # I/M, ExtP1
  c(0,1,-1,0,0,0,0,0,0,0,0,0), # M/R, ExtP3
  c(0,0,0,0,0,1,-1,0,0,0,0,0), # M/R, FlexP3
  c(0,0,0,0,0,0,0,0,0,1,-1,0), # M/R, ExtP1
  c(0,0,1,-1,0,0,0,0,0,0,0,0), # R/L, ExtP3
  c(0,0,0,0,0,0,1,-1,0,0,0,0), # R/L, FlexP3
  c(0,0,0,0,0,0,0,0,0,0,1,-1) # R/L, ExtP1
)
contrasts(group)<-t(ginv(M.finger))
fitM2.3.REML.finger <- lme(F ~ group,
                             random=list(individual=pdBlocked(list(pdIdent(~1),
                                                                   pdDiag(~location-1),
                                                                   pdIdent(~finger-1),
                                                                   pdIdent(~location:finger-1)))),
                             weights=varIdent(form=~1|location*Index),
                             correlation=corSymm(form=~1|individual/trial),
                             method="REML", control=lmeControl(msMaxIter=1000))
#***Table 16***
summary(fitM2.3.REML.finger)
#***