setwd("D:/CECILE_PL/CONFERENCES/2010_ATelier_INSERM/Phase_pratique/TD")


################### Importation of data tables #####################

# read file
bvrt<-read.table("paquid_long_570.csv",header=T,sep=";")
demo<-read.table("paquid_demo_570.csv",header=T,sep=";")

# edit data
edit(demo)

#number of subjects
tab<-table(bvrt[,1])
length(tab)

##################### BENTON analysis ########################


##### pre-analysis

# libraries required 
# lcmm for hlme function; lattice for plots & splines for spline regression
require(lcmm)
library(lattice)
library(splines)

# creation of time variable "delai"
bvrt$delai<-bvrt$age-bvrt$age0

# see the table
edit(bvrt)

# summary of the table
summary(bvrt)

# spaghetti plots 
xyplot(benton ~ delai, bvrt, groups = NUMERO, type = c('g','l'), aspect = 'xy') 

# creation of covariates variables
bvrt$delai2<-bvrt$delai*bvrt$delai/100
bvrt$cep<-1-bvrt$cep
bvrt$cepdel<-bvrt$cep*bvrt$delai
bvrt$cepdel2<-bvrt$cepdel*bvrt$delai/100
bvrt$primo<-ifelse(bvrt$age==bvrt$age0,1,0)
bvrt$age65<-bvrt$age0-65
bvrt$age65del<-bvrt$age65*bvrt$delai
bvrt$age65del2<-bvrt$age65del*bvrt$delai/100
bvrt$maledel<-bvrt$male*bvrt$delai
bvrt$maledel2<-bvrt$maledel*bvrt$delai/100


###### linear mixed model: choice of shape of trajectory
# + random effects + primo passation as single covariate
		

# no random effect: L= -6127.61
m1<-hlme(benton~delai+delai2+primo,
			random=~-1,subject=NUMERO,ng=1,data=bvrt)
m1
# random intercept: L=-5639.03
m1<-hlme(benton~delai+delai2+primo,
			random=~1, subject=NUMERO,ng=1,data=bvrt)
m1
# random intercept, delai: L= -5629.73
m1<-hlme(benton~delai+delai2+primo,
			random=~1+delai, subject=NUMERO,ng=1,data=bvrt)
m1

# random intercept, delai & delai2: L=-5625.9 
m2<-hlme(benton~delai+delai2+primo,
			random=~delai+delai2,subject=NUMERO,ng=1,data=bvrt)

m2
summary(m2)



##### Quadratic model: choice for covariate adjustment
# all: -5536.06
m2<-hlme(benton~delai+delai2+primo+cep+cepdel+cepdel2+age65+age65del+age65del2
			+male+maledel+maledel2,
			random=~delai+delai2,subject=NUMERO,ng=1,data=bvrt)
m2
summary(m2)
# without cep*t,t2: L=-5536.24
m2<-hlme(benton~delai+delai2+primo+cep+age65+age65del+age65del2
			+male+maledel+maledel2,
			random=~delai+delai2,subject=NUMERO,ng=1,data=bvrt)
# no male*t,t2 : L=-5536.96 
m2<-hlme(benton~delai+delai2+primo+cep+age65+age65del+age65del2+male,
			random=~delai+delai2,subject=NUMERO,ng=1,data=bvrt)
# no interaction with t,t2 : L=-5541.97 
m2<-hlme(benton~delai+delai2+primo+cep+age65
			+male,
			random=~delai+delai2,subject=NUMERO,ng=1,data=bvrt)
# without age65*t : L=-5537.05 
# without age65*t2: L=-5537.29 
m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
			random=~delai+delai2,subject=NUMERO,ng=1,data=bvrt)

m2
summary(m2)
plot(m2)

#### Predictions 

# table for predictions
datapred<-data.frame(intercept=rep(1,100),delai=seq(0,15,length=100)
,delai2=rep(0,100)
,primo=rep(0,100),cep=rep(0,100),male=rep(0,100),age65=rep(0,100))
		
datapred$delai2<-datapred$delai*datapred$delai/100
datapred$age65del2<-datapred$delai*datapred$age65*datapred$delai/100

# for G=1
predG1<-datapred
predG1$pred<-12.90691036-0.05679149 *predG1$delai+ 0.09168697 *predG1$delai2
plot(predG1$pred~predG1$delai,xlim=c(0,15),ylim=c(0,15),type='l',ylab="predicted benton",lwd=1,lty=1)


##### heterogeneous profiles 

### nG=1 :   BIC: 11169.29  

### nG=2
m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=2,data=bvrt)
postprob(m2)
#  BIC: 11127.78 
  

m2
summary(m2)
plot.predict.hlme(m2,datapred,"delai")

m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=2,data=bvrt
	,B=c(0,15,10,-0.2, -1, 0,0.1, -0.69531081 ,-1.75028554,
	 0.57673270,-0.09253981,-0.04735902,2.8831338,-0.1188831,
	 0.06396342,0.7405395, -0.41164459, 2.951428,1.626763))
#  BIC: 11127.78 

m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=2,data=bvrt
	,B=c(0,15,14,-0.2, -0.5, 0,0, -0.69531081 ,-1.75028554,
	 0.57673270,-0.09253981,-0.04735902,2.8831338,-0.1188831,
	 0.06396342,0.7405395, -0.41164459, 2.951428,1.626763))
m2
#   BIC: 11129.61  


m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=2,data=bvrt
	,B=c(0,10,10,-0.2, -0.5, 0,1,-0.69531081 ,-1.75028554,
	 0.57673270,-0.09253981,-0.04735902,2.8831338,-0.1188831,
	 0.06396342,0.7405395, -0.41164459, 2.951428,1.626763))
m2
#   BIC: 11129.61  


m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=2,data=bvrt
	,B=c(1,13,15,-0.1,-1,-0.1,0.1,-0.69531081 ,-1.75028554,
	 0.57673270,-0.09253981,-0.04735902,2.8831338,-0.1188831,
	 0.06396342,0.7405395, -0.41164459, 2.951428,1.626763))
m2

#  BIC: 11129.61 

m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=2,data=bvrt
	,B=c(1,13,15,-0.1,-1,-0.1,0.1,0,0,0,0,0, 2 
	,0,0.05,0,0,2,1.6))
m2
#   BIC: 11129.61  
summary(m2)
postprob(m2)


##### cubic trajectories ??
bvrt$delai3<-bvrt$delai2*bvrt$delai
# G=1  BIC: 11174.87 
m22<-hlme(benton~delai+delai2+delai3+primo+cep+male+age65+age65del2,		
	random=~delai+delai2,subject=NUMERO,ng=1,data=bvrt)
m22
# G=2  BIC: 11138.07  
m22<-hlme(benton~delai+delai2+delai3+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2+delai3,		
	random=~delai+delai2,subject=NUMERO,ng=2,data=bvrt)

# prediction avec delai 3
datapred3<-data.frame(intercept=rep(1,100),delai=seq(0,15,length=100)
,delai2=rep(0,100),delai3=rep(0,100)
,primo=rep(0,100),cep=rep(0,100),male=rep(0,100),age65=rep(0,100))
		
datapred3$delai2<-datapred3$delai*datapred3$delai/100
datapred3$delai3<-datapred3$delai2*datapred3$delai
datapred3$age65del2<-datapred3$delai*datapred3$age65*datapred3$delai/100
plot.predict.hlme(m22,datapred3,"delai")

postprob(m22)
postprob(m2)

# G=3  BIC: 11126.11 
m32<-hlme(benton~delai+delai2+delai3+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2+delai3,		
	random=~delai+delai2,subject=NUMERO,ng=3,data=bvrt)
m32
plot.predict.hlme(m32,datapred3,"delai")

postprob(m32)

# G=4  BIC: 11143.76 
m42<-hlme(benton~delai+delai2+delai3+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2+delai3,		
	random=~delai+delai2,subject=NUMERO,ng=4,data=bvrt)
m42
plot.predict.hlme(m42,datapred3,"delai")
postprob(m42)



####### NG=3 

m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=3,data=bvrt
	,B=c(1,1,13,15,10,-0.1,-1,0,-0.1,0.1,0,-0.69531081 ,-1.75028554,
	 0.57673270,-0.09253981,-0.04735902,2.8831338,-0.1188831,
	 0.06396342,0.7405395, -0.41164459, 2.951428,1.626763))

m2
summary(m2)
plot.predict.hlme(m2,datapred,"delai")
#  BIC: 11125.56 


m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=3,data=bvrt)
postprob(m2)
#   BIC: 11154.99  


m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=3,data=bvrt
	,B=c(0,0,13,15,10,-0.1,-0.5,-0.4,-1,-3,-1,-0.69531081 ,-1.75028554,
	 0.57673270,-0.09253981,-0.04735902,2.8831338,-0.1188831,
	 0.06396342,0.7405395, -0.41164459, 2.951428,1.626763))
m2
summary(m2)
#  BIC: 11125.56 

m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=3,data=bvrt
	,B=c(0,0,13,13,13,-0.1,-0.5,-0.4,-1,-2,-3,-0.69531081 ,-1.75028554,
	 0.57673270,-0.09253981,-0.04735902,2.8831338,-0.1188831,
	 0.06396342,0.7405395, -0.41164459, 2.951428,1.626763))

m2
summary(m2)
plot.predict.hlme(m2,datapred,"delai")
postprob(m2)
#     maximum log-likelihood: -5489.81  
#     AIC: 11025.61  
#     BIC: 11125.56 
summary(m2)


##### primo passation
m2<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2+primo,		
	random=~delai+delai2,subject=NUMERO,ng=3,data=bvrt
	,B=c(0,0,8,13,10,-0.1,0,0,-0.1,0.5,0.5,-0.69531081,-0.69531081,-0.69531081 ,-1.75028554,
	 0.57673270,-0.09253981,-0.04735902,2.8831338,-0.1188831,
	 0.06396342,0.7405395, -0.41164459, 2.951428,1.626763))
m2
summary(m2)
plot.predict.hlme(m2,datapred,"delai")
postprob(m2)

#     maximum log-likelihood: -5487.19  
#     AIC: 11024.38  
#     BIC: 11133.02
 
prob<-m2$pprob
prob2<-cbind(prob,demo)

summary(prob2)

## significant association: demented,cep,mvie
## negative association: male,age65,SUJDEPR0,SUJPSY0,veuf0,gseul0,fdrcv0
crois<-table(prob2$class,prob2$sommed0)
crois
chisq.test(crois)


#### NG=4
m3<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=4,data=bvrt
	,B=c(0,0,0,8,13,10,14,-0.1,-0.2,0.3,-0.4,-1,-2,-3,-4,-0.69531081,-1.75028554,
	 0.57673270,-0.09253981,-0.04735902,2.8831338,-0.1188831,
	 0.06396342,0.7405395, -0.41164459, 2.951428,1.626763))

m3
summary(m3)
postprob(m3)
plot.predict.hlme(m3,datapred,"delai")
#   BIC: 11138.16  

m3<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=4,data=bvrt)
m3
#  BIC: 11138.28  


m3<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	mixture=~delai+delai2,		
	random=~delai+delai2,subject=NUMERO,ng=4,data=bvrt
	,B=c(0,0,0,8,10,13,14,-0.4,-0.1,-0.2,-0.3,-2,-4,-1,-3,-0.69531081,-1.75028554,
	 0.57673270,-0.09253981,-0.04735902,2.8831338,-0.1188831,
	 0.06396342,0.7405395, -0.41164459, 2.951428,1.626763))

m3
summary(m3)
#   BIC: 11144.88  




###################### GROUP-BASED TRAJECTORY #####################

GBT1<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,	
	random=~-1,
	subject=NUMERO,ng=1,data=bvrt)
GBT1
# G=1   maximum log-likelihood: -5879.74  
#     AIC: 11777.48  
#     BIC: 11816.59  
GBT1<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	random=~-1,
	mixture=~delai+delai2,		
	subject=NUMERO,ng=2,data=bvrt)
GBT1
# G=2       maximum log-likelihood: -5617.69  
#     AIC: 11261.38  
#     BIC: 11317.88  
GBT1<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	random=~-1,
	mixture=~delai+delai2,		
	subject=NUMERO,ng=3,data=bvrt)
GBT1
summary(GBT1)
# G=3     maximum log-likelihood: -5531.22  
     AIC: 11096.44  
     BIC: 11170.31 
GBT1<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	random=~-1,
	mixture=~delai+delai2,		
	subject=NUMERO,ng=4,data=bvrt)
GBT1
# G=4         maximum log-likelihood: -5502.48  
     AIC: 11046.97  
     BIC: 11138.22  
GBT1<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	random=~-1,
	mixture=~delai+delai2,		
	subject=NUMERO,ng=5,data=bvrt)
GBT1
# G=5         maximum log-likelihood: -5490.91  
     AIC: 11031.82  
     BIC: 11140.46  
GBT1<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	random=~-1,
	mixture=~delai+delai2,		
	subject=NUMERO,ng=6,data=bvrt)
GBT1
# G=6  maximum log-likelihood: -5484.31  
     AIC: 11026.63  
     BIC: 11152.65  
GBT1<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	random=~-1,
	mixture=~delai+delai2,		
	subject=NUMERO,ng=6,data=bvrt,
	B=c(0,0,0,0,0,15,13,13,10,10,7,0,-0.5,0,-0.2,-0.3,0,-0.2,
		-0.4,0,-0.5,-1,-1,-0.7,-1.2,0,-0.2,-0.2,2))
GBT1
# G=6  maximum log-likelihood: -5484.31  
     AIC: 11026.63  
     BIC: 11152.65  

GBT1<-hlme(benton~delai+delai2+primo+cep+male+age65+age65del2,
	random=~-1,
	mixture=~delai+delai2,		
	subject=NUMERO,ng=6,data=bvrt,
	B=c(0,0,0,0,0,15,13,13,10,10,7,0,-0.5,0,-0.2,0.3,0,0.2,
		-0.4,0,-0.5,-1,-1,-0.7,-1.2,0,-0.2,-0.2,2))
GBT1