No need code of dementia.Rmd

---
title: "No need code"
output: html_document
date: "2023-06-29"
---


```{r (need to change) only keep those have mmse follow-up}

data.nommse.follow<-data1%>%
  filter(is.na(mmseyr1total)&is.na(mmseyr2total)&is.na(mmseyr3total)&is.na(mmseyr4total)&is.na(mmseyr5total)&is.na(mmseyr6total)&is.na(mmseyr7total)&is.na(mmseyr8total)&is.na(mmseyr9total)&is.na(mmseyr10total)&is.na(mmseyr11total)&is.na(mmseyr12total))

data.nommse.follow%>%
  group_by(study)%>%
  count()

#ICICLE 14, NYPUM 6, ParkWest 3, PICNICS 19, PINE 5 (P158,P529,P579,P794,P85)
#47 patients didn't have mmse follow-up 

data1<-data1%>%
  filter(!is.na(mmseyr1total)|!is.na(mmseyr2total)|!is.na(mmseyr3total)|!is.na(mmseyr4total)|!is.na(mmseyr5total)|!is.na(mmseyr6total)|!is.na(mmseyr7total)|!is.na(mmseyr8total)|!is.na(mmseyr9total)|!is.na(mmseyr10total)|!is.na(mmseyr11total)|!is.na(mmseyr12total))

6+3+19+5 #Without ICICIE 

data1%>%
  group_by(study)%>%
  count()

119+130+185+230+275 #Without ICICLE 939
```

```{r not run match longtidue data with baseline1}

#If patients only had one follow-up, which mean the follow-up equal to the event year, then need to be removed. 

#setdiff(unique(data.baseline.1$id),unique(data_long.1$id)) #drop off the year0 then can see #65 182 351 689 737 839

# List of dataset names
dataset_names <- c("data.baseline.1", "data.baseline.2", "data.baseline.3", "data.baseline.4", "data.baseline.5", "data.baseline.6", "data.baseline.7", "data.baseline.8","data.baseline.9")

# List of ids to be removed
ids_to_remove <- c(65, 182, 351, 689, 737, 839)

# Loop through each dataset
for (i in seq_along(dataset_names)) {
  # Get the dataset
  dataset <- get(dataset_names[i])
  
  # Remove specific ids from the dataset
  for (id in ids_to_remove) {
    dataset <- dataset[dataset$id != id, ]
  }
  
  # Assign the result to the output dataset
  assign(dataset_names[i], dataset)
}


data.baseline.1%>%
  group_by(study)%>%
  count()

```



```{r Not run event years should be longer than followupyears}

data_long.1<-data_long.1%>%
  filter(years>followupyears)

data_long.2<-data_long.2%>%
  filter(years>followupyears)

data_long.3<-data_long.3%>%
  filter(years>followupyears)

data_long.4<-data_long.4%>%
  filter(years>followupyears)

data_long.5<-data_long.5%>%
  filter(years>followupyears)

data_long.6<-data_long.6%>%
  filter(years>followupyears)

data_long.7<-data_long.7%>%
  filter(years>followupyears)

data_long.8<-data_long.8%>%
  filter(years>followupyears)

data_long.9<-data_long.9%>%
  filter(years>followupyears)

setdiff(unique(data.baseline.1$idpicc),unique(data_long.1$idpicc)) #internal cens? Got dementia between baseline the next follow-up time

data.baseline.1%>%
  filter(idpicc %in% c(2,6,7,15,22,23,32,34,45,50,59,62,65,66,78,92,94,99,101,314,343,405,456,
                       634,637,644,646,653,661,662,664,665,666,676,680,683,686,689,690,698,705,712,716,717,719,721,
                       734,735,756,759,760,773,781,783,785,787,792,801,804,807,818,831,872,874,895,896,903,909,1577,
                       1146,1599,1155,1158,1171,1639,1645,1666,1677,1679,1313,1687,1588))%>%
  group_by(study)%>%
  count()


# List of dataset names
dataset_names <- c("data.baseline.1", "data.baseline.2", "data.baseline.3", "data.baseline.4", "data.baseline.5", "data.baseline.6", "data.baseline.7", "data.baseline.8","data.baseline.9")

# List of ids to be removed
idpiccs_to_remove <- c(2,6,7,15,22,23,32,34,45,50,59,62,65,66,78,92,94,99,101,314,343,405,456,
                       634,637,644,646,653,661,662,664,665,666,676,680,683,686,689,690,698,705,712,716,717,719,721,
                       734,735,756,759,760,773,781,783,785,787,792,801,804,807,818,831,872,874,895,896,903,909,1577,
                       1146,1599,1155,1158,1171,1639,1645,1666,1677,1679,1313,1687,1588)

# Loop through each dataset
for (i in seq_along(dataset_names)) {
  # Get the dataset
  dataset <- get(dataset_names[i])
  
  # Remove specific ids from the dataset
  for (idpicc in idpiccs_to_remove) {
    dataset <- dataset[dataset$idpicc != idpicc, ]
  }
  
  # Assign the result to the output dataset
  assign(dataset_names[i], dataset)
}


data.baseline.1%>%
  group_by(study)%>%
  count()

data.baseline.1%>%
  count()

```



```{r Characteristic baseline}


# the id need to remove before c(65, 182, 351, 689, 737, 839)

imp3.new[imp3.new$id==c(65, 182, 351, 689, 737, 839),]%>%
  select(study,idpicc,originalid)

#study      idpicc originalid
#CamPalGN    69    244
#NYPUM      366    337
#ParkWest   546    B100
#PINE       1600   P193
#           1188   P351
#           1271   P659

data1.10.remove<-data1.10[data1.10$idpicc!=69,]
data1.10.remove<-data1.10.remove[data1.10.remove$idpicc!=366,]
data1.10.remove<-data1.10.remove[data1.10.remove$idpicc!=546,]
data1.10.remove<-data1.10.remove[data1.10.remove$idpicc!=1600,]
data1.10.remove<-data1.10.remove[data1.10.remove$idpicc!=1188,]
data1.10.remove<-data1.10.remove[data1.10.remove$idpicc!=1271,]

data1.10.remove<-data1.10.remove%>%
  filter(! idpicc %in% c(2,6,7,15,22,23,32,34,45,50,59,62,65,66,78,92,94,99,101,314,343,405,456,
                       634,637,644,646,653,661,662,664,665,666,676,680,683,686,689,690,698,705,712,716,717,719,721,
                       734,735,756,759,760,773,781,783,785,787,792,801,804,807,818,831,872,874,895,896,903,909,1577,
                       1146,1599,1155,1158,1171,1639,1645,1666,1677,1679,1313,1687,1588))
```


```{r missing pattern}

data.baseline.noIC<-data.baseline.new

data.baseline.noIC%>%
  group_by(study)%>%
  count()


data.rename.noIC<-data.baseline.noIC%>%
  rename("Age at baseline"=agebl, 
         "Sex"=sex,
         "MDS-UPDRS part3"=mdsupdrspart3bltotalconvertedasa,"Hoehn and Yahr Scale"=hybl, "Smoking"=smoking,
         "Hallucinations"=hallucinationsindex,"Cognitive symptoms"= cognitiveindex
         ) #

explanatory<-c("Smoking","MDS-UPDRS part3","Hoehn and Yahr Scale","Hallucinations","Cognitive symptoms") 
dependent<- c("cens","tt") #"Smoking"

mispattern<-data.rename.noIC %>% 
  missing_pattern(explanatory)

#---check how many missing at least one---

colnames(data.baseline.noIC)

nrow(data.baseline.noIC[complete.cases(data.baseline.noIC[,7:10]),]) 

896-825
71/896

```

```{r jomo exclude ICICLE}

#-----EXCLUDE ICICLE----


data.baseline.noIC$cons<-1

data.baseline.noIC$nelsonaalen<-nelsonaalen(data.baseline.noIC,years,cens) #0 is right censor,1 is event

Y<- data.baseline.noIC[,c("mdsupdrspart3bltotalconvertedasa","hallucinationsindex","cognitiveindex","smoking")]  


X<-data.baseline.noIC[,c("cons","agebl","sex","hybl","nelsonaalen")] #adding Nelson-Aalen estimate  

clus<-data.baseline.noIC$study

imp.dry<-jomo.MCMCchain(Y = Y,X = X,clus = clus, nburn = 2)

set.seed(15678)
imp1 <- jomo.MCMCchain(Y = Y, X = X, clus = clus, nburn = 5000)


#head(imp1$collectbeta) # check beta


#plot trace for each parameter value



#png("Jomo1.png",width = 3500,height =2000,res = 400)

#par(mfrow=c(2,4))

plot(imp1$collectbeta[1, 1, 1:5000], type = "l", ylab = expression(beta["mdsupdrspart3,0"]),
     xlab = "Iteration number" )

plot(imp1$collectbeta[1, 2, 1:5000], type = "l", ylab = expression(beta["hallucinationsindex.1 ,0"]),
     xlab = "Iteration number" )

plot(imp1$collectbeta[1, 3, 1:5000], type = "l", ylab = expression(beta["cognitiveindex.1,0"]),
     xlab = "Iteration number" ) 

plot(imp1$collectbeta[1, 4, 1:5000], type = "l", ylab = expression(beta["cognitiveindex.2,0"]),
     xlab = "Iteration number" ) 

plot(imp1$collectbeta[1, 5, 1:5000], type = "l", ylab = expression(beta["smoking.1,0"]),
     xlab = "Iteration number" ) 

plot(imp1$collectbeta[1, 6, 1:5000], type = "l", ylab = expression(beta["smoking.2,0"]),
     xlab = "Iteration number" ) 

#plot trace for cov matrix element
#imp1$collectomega[,,1] #check the row and col name
#Category variable don't need to plot, just a straight line

plot(imp1$collectomega[1, 1, 1:5000], type = "l", ylab = expression(omega[MDS-UPDRS,1,1]^2),
     xlab = "Iteration number" )


#dev.off()


# Capture the state of the sampler as starting values for the second set of iterations:
beta.start <- imp1$collectbeta[,,5000] # capture the fixed parameter values
l1cov.start <- imp1$collectomega[,,5000] # capture the level-1 covariance matrix values
start.imp <- imp1$finimp.latnorm # capture the final imputed data set 



#Re-run the same function for a larger number of iterations
imp2 <- jomo.MCMCchain(Y = Y, X = X, clus = clus, beta.start = beta.start, l1cov.start = l1cov.start,
                       start.imp = start.imp, nburn = 5000)

# Check the trace again

#png("Jomo2.png",width = 3500,height =2000,res = 400)

#par(mfrow=c(2,4))

plot(imp2$collectbeta[1, 1, 1:5000], type = "l", ylab = expression(beta["mdsupdrspart3,0"]),
     xlab = "Iteration number" )

plot(imp2$collectbeta[1, 2, 1:5000], type = "l", ylab = expression(beta["hallucinationsindex.1,0"]),
     xlab = "Iteration number" ) 

plot(imp2$collectbeta[1, 3, 1:5000], type = "l", ylab = expression(beta["cognitiveindex.1,0"]),
     xlab = "Iteration number" ) 

plot(imp2$collectbeta[1, 4, 1:5000], type = "l", ylab = expression(beta["cognitiveindex.2,0"]),
     xlab = "Iteration number" ) 

plot(imp2$collectbeta[1, 5, 1:5000], type = "l", ylab = expression(beta["smoking.1,0"]),
     xlab = "Iteration number" ) 

plot(imp2$collectbeta[1, 6, 1:5000], type = "l", ylab = expression(beta["smoking.2,0"]),
     xlab = "Iteration number" ) 

#plot trace for cov matrix element

plot(imp2$collectomega[1, 1, 1:5000], type = "l", ylab = expression(omega[mdsupdrspart3,1,1]^2),
     xlab = "Iteration number" )



#dev.off()

#collect posterior mean of cov matrix
l1cov.guess <- apply(imp2$collectomega, c(1, 2), mean)

dim(imp2$collectomega[,,1]) #degrees of freedom is 6

# Multiply by degrees of freedom to get scale matrix

l1cov.prior <- l1cov.guess*6

# Perform multilevel imputation:
imp3 <- jomo(Y = Y, X = X, clus = clus, l1cov.prior = l1cov.prior, nburn = 5000, nbetween = 1000, nimp = 7,meth = "random" )  #The nimp is the precentage of missing 


```


```{r merge data baseline}

names(imp3)
dim(imp3)
#View original (partially observed) data:
head(imp3)
# View last imputation (the left most column is the row number):
head(imp3[imp3$Imputation == 1,])

#create id to merge the data 

data.baseline$id<-seq(nrow(data.baseline)) 

data.baseline.time<-data.baseline%>%
  select(id,study,idpicc,cens,tt,years)

imp3.new<-merge(imp3,data.baseline.time,by.x = "id",by.y = "id")

```


```{r check the imputation set}

imp3.new$age10<-imp3.new$agebl/10

imp3.new$mdsupdrs3.10<-imp3.new$mdsupdrspart3bltotalconvertedasa/10


imp3.new[imp3.new$mdsupdrs3.10<0,]

#Just one imputated value is negative

imp3.new[imp3.new$id==435&imp3.new$Imputation==6,]$mdsupdrs3.10<-0

summary(imp3.new)  

data.baseline.imp<-imputationList(split(imp3.new, imp3.new$Imputation)[-1])
```


```{r just one study}


data.lm.cam<-data.lm.new%>%
  filter(study=="CamPalGN")

data.baseline.cam<-data.baseline.imp[["imputations"]][["1"]]%>%
  filter(study=="CamPalGN")


length(unique(data.baseline.cam$idpicc)) #119
length(unique(data.lm.cam$idpicc)) #118

unique(data.baseline.cam$idpicc)
unique(data.lm.cam$idpicc)
#data.baseline.cam has 69, but data.lm.cam doesn't

data1%>%
  filter(idpicc==69)

data_long%>%
  filter(idpicc==69)

#This is because the followupyears and years are the same in 69 patients. Patients death soon after the only time followup view at year 3.

data.baseline.cam.new<-subset(data.baseline.cam,unique(data.baseline.cam$idpicc) %in% unique(data.lm.cam$idpicc))

baseline.cox.cam<-coxph(Surv(years,cens)~age10+sex+mdsupdrs3.10+hybl+hallucinationsindex+cognitiveindex,data=data.baseline.cam.new, x = TRUE,model = TRUE)

lmeFit.cam<- lme(mmse~ followupyears, random = ~ followupyears | idpicc, data = data.lm.cam)

jm(baseline.cox.cam, lmeFit.cam, time_var = "followupyears")
```


```{r baseline Cox regression chose one to run}



#imp1.data<-data.baseline.imp[["imputations"]][["1"]]

#imp1.data<-subset(imp1.data,unique(data.baseline.noIC$idpicc) %in% unique(data.lm$idpicc))

baseline.cox1<-coxph(Surv(years,cens)~strata(clus)+age10+sex+mdsupdrs3.10+hybl+hallucinationsindex+cognitiveindex,data=data.baseline.imp[["imputations"]][["1"]],cluster=idpicc, x = TRUE,model = TRUE)

```


```{r try randome effect linear}

ggplot(data = mmse_plot, aes(x = followupyears, y = mmse))+
    geom_point()+
    geom_smooth(method = "loess")

#---exclude ICICLE


data.lm<-data.lm%>%
  rename("clus"=study)


data.lm%>%
  filter(is.na(mmse))

#remove missing mmse in the data

#data.lm<-data.lm[complete.cases(data.lm),]


lmeFit<- lme(mmse~ followupyears+followupyears:clus, random = ~ followupyears | idpicc, data = data.lm)


```


```{r without loop}




data.imp1$year1<-as.numeric(as.Date(as.character(data.imp1$datevisityr1), format="%Y-%m-%d")-
                as.Date(as.character(data.imp1$datevisitbl), format="%Y-%m-%d"))/365.25

data.imp1$year2<-as.numeric(as.Date(as.character(data1.10$datevisityr2), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25


data1.10$year3<-as.numeric(as.Date(as.character(data1.10$datevisityr3), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year4<-as.numeric(as.Date(as.character(data1.10$datevisityr4), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year5<-as.numeric(as.Date(as.character(data1.10$datevisityr5), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year6<-as.numeric(as.Date(as.character(data1.10$datevisityr6), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year7<-as.numeric(as.Date(as.character(data1.10$datevisityr7), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year8<-as.numeric(as.Date(as.character(data1.10$datevisityr8), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year9<-as.numeric(as.Date(as.character(data1.10$datevisityr9), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year10<-as.numeric(as.Date(as.character(data1.10$datevisityr10), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year11<-as.numeric(as.Date(as.character(data1.10$datevisityr11), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year12<-as.numeric(as.Date(as.character(data1.10$datevisityr12), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year13<-as.numeric(as.Date(as.character(data1.10$datevisityr13), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year14<-as.numeric(as.Date(as.character(data1.10$datevisityr14), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25

data1.10$year15<-as.numeric(as.Date(as.character(data1.10$datevisityr15), format="%Y-%m-%d")-
                as.Date(as.character(data1.10$datevisitbl), format="%Y-%m-%d"))/365.25


```



```{r Random intercept model}

#Random intercept model (each study have random intercept)

class(data_long.1$study)

Mlme1.1<-lme(mmse~ followupyears, random = ~1 | study, data=data_long.1)

summary(Mlme1.1)

F0<-fitted(Mlme1.1,level = 0) #fitted values obtained by the population model
F1<-fitted(Mlme1.1,level = 1) #gives the within-study values


study_colors <- c("blue", "red", "green", "purple", "pink")
study_names<-c("CamPalGN","NYPUM","ParkWest","PICNICS","PINE")
study_color_mapping <- setNames(study_colors, study_names)

plot(sort(data_long.1$followupyears), F0 [order(data_long.1$followupyears)],lwd=2,type = "l",ylab = "mmse",xlab = "follow-up years")
for (i in study_names) {
  x1<-data_long.1 [data_long.1$study==i,]$followupyears
  y1<-F1[data_long.1$study == i]
  lines(sort(x1),y1[order(x1)],col=study_color_mapping[i],lwd=2)
}

legend("topright", legend = study_names, col=study_colors,lwd=2)



#Random intercept model (each patients have random intercept)

Mlme2.1<-lme(mmse~ followupyears, random = ~1 | idpicc, data=data_long.1)

summary(Mlme2.1) #21349.39 21374.4


#random intercept for study

Mlme1.1<-lme(mmse~ followupyears, random = ~1 | study, data=data_long.1)

#random slope for 

Mlme1.1<-lme(mmse~ followupyears, random = ~1 + followupyears | study, data=data_long.1)


#random intercept for patients


Mlme3.1<-lme(mmse~ followupyears, random = ~ 1 | idpicc, data=data_long.1)

summary(Mlme1.1) #random intecept for study
summary(Mlme3.1) #random intercept for patients
summary(Mlme4.1) #random intercept for patients and random slope for followupyears
summary(Mlme5.1) #random intercept for patients




#random slope for followupyears and random slope for patients

Mlme4.1<-lme(mmse~ followupyears , random = ~ followupyears | idpicc, data=data_long.1)
summary(Mlme4.1) #19644.03 19681.55

#mmse outcome, study as factor, followupyears


Mlme5.1<-lme(mmse~ followupyears +study , random = ~ followupyears | idpicc, data=data_long.1)
summary(Mlme5.1)
```


```{r Not run combine baseline dataset again to run Cox}

combined_data.baseline<-rbind(data.baseline.1,data.baseline.2,data.baseline.3,data.baseline.4,data.baseline.5,data.baseline.6,data.baseline.7,data.baseline.8,data.baseline.9)

imp.baseline.list<-jomo2mitml.list(combined_data.baseline)

try.cox<-with(data=imp.baseline.list,coxph(Surv(years,cens)~age10+sex+yearseducation+mdsupdrs3.10+hybl+hallucinationsindex+cognitiveindex+strata(study)),x=T,y=T)

```

```{r Not run combine long data to run lme}

combined_data_long<-rbind(data_long.1,data_long.2,data_long.3,data_long.4,data_long.5,data_long.6,data_long.7,data_long.8,data_long.9)

imp.data_long<-jomo2mitml.list(combined_data_long)

try.lme<-with(data=imp.data_long,lme(mmse~ agebl+yearseducation+followupyears, random = ~ followupyears | idpicc, data=combined_data_long))

#I don't know how to put these two together.
```


```{r }

try.lme<-lme(mmse~ agebl+yearseducation+followupyears, random = ~ followupyears | idpicc,data=combined_data_long) 

try.cox<-coxph(Surv(years,cens)~age10+sex+yearseducation+mdsupdrs3.10+hybl+hallucinationsindex+cognitiveindex+strata(study),data=combined_data.baseline, x = TRUE,model = TRUE)

try.cox<-with(data=imp.baseline.list,coxph(Surv(years,cens)~age10+sex+yearseducation+mdsupdrs3.10+hybl+hallucinationsindex+cognitiveindex+strata(study)),x=T,y=T)

try.cox.res<-testEstimates(try.cox)

try.cox.estimate<-as.data.frame(try.cox.res[["estimates"]])

try.cox.estimate$Estimate



pool(try.cox)

try.jm<-jm(try.cox, try.lme, time_var = "followupyears")



#Error in checkForRemoteErrors(val) : 
#3 nodes produced errors; first error: addition: incompatible matrix dimensions: 890x1 and 8010x1

#The problem is in COX, The N=8010 (890*9). I need to figure out how to pool an average model from it
```