rm(list=ls())
#########Load packages
library(robustlmm)
library(mice)

#Load data
A<-read.csv("C:/Users/USER/Desktop/missing.csv",header=TRUE)
names(A)
head(A)
tail(A)

A$enterdate1[A$enterdate1==" "]<- -9
A$enterdate2[A$enterdate2==" "]<- -9
A$enterdate3[A$enterdate3==" "]<- -9
A$enterdate4[A$enterdate4==" "]<- -9
A$enterdate5[A$enterdate5==" "]<- -9
A$enterdate6[A$enterdate6==" "]<- -9

####################################################
#first inclusion criterion
####################################################

n<-length(A$id);n
a<-rep(NA,nrow(A))

for (i in 1:n) {
  if (A$age1[i]!=-9 | A$sex1[i]!=-9 | A$enterdate1[i]!=-9) {
a[i]=1
}
  if (A$age1[i]==-9 & A$sex1[i] ==-9 & A$enterdate1[i]==-9) {
    if (A$age2[i]!=-9| A$sex2[i]!=-9 | A$enterdate2[i]!=-9) {

a[i]=2
}
  }

  if (A$age1[i]==-9 & A$sex1[i] ==-9 & A$enterdate1[i]==-9 &A$age2[i]==-9 & A$sex2[i] ==-9 & A$enterdate2[i]==-9) {
    if (A$age3[i]!=-9| A$sex3[i]!=-9 | A$enterdate3[i]!=-9) {

a[i]=3
}
  }

  if (A$age1[i]==-9 & A$sex1[i] ==-9 & A$enterdate1[i]==-9 & A$age2[i]==-9 & A$sex2[i] ==-9 & A$enterdate2[i]==-9& A$age3[i]==-9 & A$sex3[i] ==-9 & A$enterdate3[i]==-9) {
    if (A$age4[i]!=-9| A$sex4[i]!=-9 | A$enterdate4[i]!=-9) {

a[i]=4
}
  }

  if (A$age1[i]==-9 & A$sex1[i] ==-9 & A$enterdate1[i]==-9 & A$age2[i]==-9 & A$sex2[i] ==-9 & A$enterdate2[i]==-9 & 
A$age3[i]==-9 & A$sex3[i] ==-9 & A$enterdate3[i]==-9 & A$age4[i]==-9 & A$sex4[i] ==-9 & A$enterdate4[i]==-9) {
    if (A$age5[i]!=-9| A$sex5[i]!=-9 | A$enterdate5[i]!=-9) {

a[i]=5
}
  }

  if (A$age1[i]==-9 & A$sex1[i] ==-9 & A$enterdate1[i]==-9 & A$age2[i]==-9 & A$sex2[i] ==-9 & A$enterdate2[i]==-9 & 
A$age3[i]==-9 & A$sex3[i] ==-9 & A$enterdate3[i]==-9 & A$age4[i]==-9 & A$sex4[i] ==-9 & A$enterdate4[i]==-9
& A$age5[i]==-9 & A$sex5[i] ==-9 & A$enterdate5[i]==-9) {
    if (A$age6[i]!=-9| A$sex6[i]!=-9 | A$enterdate6[i]!=-9) {

a[i]=6
}
  }
 if (A$age1[i]!=-9 & A$sex1[i] !=-9 & A$enterdate1[i]!=-9 & A$age2[i]!=-9 & A$sex2[i] !=-9 & A$enterdate2[i]!=-9 & 
A$age3[i]!=-9 & A$sex3[i] !=-9 & A$enterdate3[i]!=-9 & A$age4[i]!=-9 & A$sex4[i] !=-9 & A$enterdate4[i]!=-9
& A$age5[i]!=-9 & A$sex5[i] !=-9 & A$enterdate5[i]!=-9 & A$age6[i]!=-9 & A$sex6[i] !=-9 & A$enterdate6[i]!=-9) {
a[i]=7
}
else if (A$age1[i]==-9 & A$sex1[i] ==-9 & A$enterdate1[i]==-9 & A$age2[i]==-9 & A$sex2[i] ==-9 & A$enterdate2[i]==-9 & 
A$age3[i]==-9 & A$sex3[i] ==-9 & A$enterdate3[i]==-9 & A$age4[i]==-9 & A$sex4[i] ==-9 & A$enterdate4[i]==-9
& A$age5[i]==-9 & A$sex5[i] ==-9 & A$enterdate5[i]==-9 & A$age6[i]==-9 & A$sex6[i] ==-9 & A$enterdate6[i]==-9) {
a[i]=0
}
}

table(a)
A<-cbind(A,a)
A<-A[which(A$a>=1),]
n<-length(A$id);n

##############################################
#seconde inclusion criterion
##############################################

b1<-rep(NA,nrow(A))
for (i in 1:n){
if(A$dbp1[i] ==-9) {
b1[i]=0
} 
else 
{b1[i]=1
}
}
table(b1)

b2<-rep(NA,nrow(A))
for (i in 1:n){
if(A$dbp2[i]==-9) {
b2[i]=0
} 
else 
{b2[i]=1
}
}
table(b2)

b3<-rep(NA,nrow(A))
for (i in 1:n){
if(A$dbp3[i]==-9) {
b3[i]=0
} 
else 
{b3[i]=1
}
}
table(b3)

b4<-rep(NA,nrow(A))
for (i in 1:n){
if(A$dbp4[i]==-9) {
b4[i]=0
} 
else 
{b4[i]=1
}
}
table(b4)

b5<-rep(NA,nrow(A))
for (i in 1:n){
if(A$dbp5[i]==-9) {
b5[i]=0
} 
else 
{b5[i]=1
}
}
table(b5)

b6<-rep(NA,nrow(A))
for (i in 1:n){
if(A$dbp6[i]==-9) {
b6[i]=0
} 
else 
{b6[i]=1
}
}
table(b6)
c<-b1+b2+b3+b4+b5+b6
A<-cbind(A,c)
A<-A[which(A$c>=2),]
n<-length(A$id);n

###################################
#Transform drug to binary variable
###################################

#drug1
for (i in 1:n){
if (A$drug1[i]==9) {
A$drug1[i]=0
} 
else if (A$drug1[i]==-9) {
A$drug1[i]=0
} 
else if (A$drug1[i]==2) {
A$drug1[i]=0
} 
else {
A$drug1[i]=1
}
}
drug1<-data.frame(A$drug1)

#drug2
for (i in 1:n){
if (A$drug2[i]==9) {
A$drug2[i]=0
} 
else if (A$drug2[i]==-9) {
A$drug2[i]=0
} 
else if (A$drug2[i]==2) {
A$drug2[i]=0
} 
else {
A$drug2[i]=1
}
}
drug2<-data.frame(A$drug2)

#drug3
for (i in 1:n){
if (A$drug3[i]==9) {
A$drug3[i]=0
} 
else if (A$drug3[i]==-9) {
A$drug3[i]=0
} 
else if (A$drug3[i]==2) {
A$drug3[i]=0
} 
else {
A$drug3[i]=1
}
}
drug3<-data.frame(A$drug3)

#drug4
for (i in 1:n){
if (A$drug4[i]==9) {
A$drug4[i]=0
} 
else if (A$drug4[i]==-9) {
A$drug4[i]=0
} 
else if (A$drug4[i]==2) {
A$drug4[i]=0
} 
else {
A$drug4[i]=1
}
}
drug4<-data.frame(A$drug4)

#drug5
for (i in 1:n){
if (A$drug5[i]==9) {
A$drug5[i]=0
} 
else if (A$drug5[i]==-9) {
A$drug5[i]=0
} 
else if (A$drug5[i]==2) {
A$drug5[i]=0
} 
else {
A$drug5[i]=1
}
}
drug5<-data.frame(A$drug5)

#drug6
for (i in 1:n){
if (A$drug6[i]==9) {
A$drug6[i]=0
} 
else if (A$drug6[i]==-9) {
A$drug6[i]=0
} 
else if (A$drug6[i]==2) {
A$drug6[i]=0
} 
else {
A$drug6[i]=1
}
}
drug6<-data.frame(A$drug6)

######################################
#dbp
######################################

# dbp1=+-10
for (i in 1:n){
if (A$drug1[i]==1&A$dbp1[i]!=-9){
A$dbp1[i]=A$dbp1[i]+10
}
else if (A$drug1[i]==0&A$dbp1[i]!=-9){
A$dbp1[i]=A$dbp1[i]-10
}
else{
A$dbp1[i]=A$dbp1[i]
}
}
dbp1<-data.frame(A$dbp1)

# dbp2=+-10
for (i in 1:n){
if (A$drug2[i]==1&A$dbp2[i]!=-9){
A$dbp2[i]=A$dbp2[i]+10
}
else if (A$drug2[i]==0&A$dbp2[i]!=-9){
A$dbp2[i]=A$dbp2[i]-10
}
else{
A$dbp2[i]=A$dbp2[i]
}
}
dbp2<-data.frame(A$dbp2)

# dbp3=+-10
for (i in 1:n){
if (A$drug3[i]==1&A$dbp3[i]!=-9){
A$dbp3[i]=A$dbp3[i]+10
}
else if (A$drug3[i]==0&A$dbp3[i]!=-9){
A$dbp3[i]=A$dbp3[i]-10
}
else{
A$dbp3[i]=A$dbp3[i]
}
}
dbp3<-data.frame(A$dbp3)

# dbp4=+-10
for (i in 1:n){
if (A$drug4[i]==1&A$dbp4[i]!=-9){
A$dbp4[i]=A$dbp4[i]+10
}
else if (A$drug4[i]==0&A$dbp4[i]!=-9){
A$dbp4[i]=A$dbp4[i]-10
}
else{
A$dbp4[i]=A$dbp4[i]
}
}
dbp4<-data.frame(A$dbp4)

# dbp5=+-10
for (i in 1:n){
if (A$drug5[i]==1&A$dbp5[i]!=-9){
A$dbp5[i]=A$dbp5[i]+10
}
else if (A$drug5[i]==0&A$dbp5[i]!=-9){
A$dbp5[i]=A$dbp5[i]-10
}
else{
A$dbp5[i]=A$dbp5[i]
}
}
dbp5<-data.frame(A$dbp5)

# dbp6=+-10
for (i in 1:n){
if (A$drug6[i]==1&A$dbp6[i]!=-9){
A$dbp6[i]=A$dbp6[i]+10
}
else if (A$drug6[i]==0&A$dbp6[i]!=-9){
A$dbp6[i]=A$dbp6[i]-10
}
else{
A$dbp6[i]=A$dbp6[i]
}
}
dbp6<-data.frame(A$dbp6)

######################################
#sbp
######################################

#sbp1=+-10
for (i in 1:n){
if (A$drug1[i]==1&A$sbp1[i]!=-9){
A$sbp1[i]=A$sbp1[i]+10
}
else if (A$drug1[i]==0&A$sbp1[i]!=-9){
A$sbp1[i]=A$sbp1[i]-10
}
else{
A$sbp1[i]=A$sbp1[i]
}
}
sbp1<-data.frame(A$sbp1)

#sbp2=+-10
for (i in 1:n){
if (A$drug2[i]==1&A$sbp2[i]!=-9){
A$sbp2[i]=A$sbp2[i]+10
}
else if (A$drug2[i]==0&A$sbp2[i]!=-9){
A$sbp2[i]=A$sbp2[i]-10
}
else{
A$sbp2[i]=A$sbp2[i]
}
}
sbp2<-data.frame(A$sbp2)

#sbp3=+-10
for (i in 1:n){
if (A$drug3[i]==1&A$sbp3[i]!=-9){
A$sbp3[i]=A$sbp3[i]+10
}
else if (A$drug3[i]==0&A$sbp3[i]!=-9){
A$sbp3[i]=A$sbp3[i]-10
}
else{
A$sbp3[i]=A$sbp3[i]
}
}
sbp3<-data.frame(A$sbp3)

#sbp4=+-10
for (i in 1:n){
if (A$drug4[i]==1&A$sbp4[i]!=-9){
A$sbp4[i]=A$sbp4[i]+10
}
else if (A$drug4[i]==0&A$sbp4[i]!=-9){
A$sbp4[i]=A$sbp4[i]-10
}
else{
A$sbp4[i]=A$sbp4[i]
}
}
sbp4<-data.frame(A$sbp4)

#sbp5=+-10
for (i in 1:n){
if (A$drug5[i]==1&A$sbp5[i]!=-9){
A$sbp5[i]=A$sbp5[i]+10
}
else if (A$drug5[i]==0&A$sbp5[i]!=-9){
A$sbp5[i]=A$sbp5[i]-10
}
else{
A$sbp5[i]=A$sbp5[i]
}
}
sbp5<-data.frame(A$sbp5)

#sbp6=+-10
for (i in 1:n){
if (A$drug6[i]==1&A$sbp6[i]!=-9){
A$sbp6[i]=A$sbp6[i]+10
}
else if (A$drug6[i]==0&A$sbp6[i]!=-9){
A$sbp6[i]=A$sbp6[i]-10
}
else{
A$sbp6[i]=A$sbp6[i]
}
}
sbp6<-data.frame(A$sbp6)

#######################################
# Data preparation
#######################################
A[A==-9]<-NA
names(A)
A<-A[,-21:-58]
names(A)
#reshape to long
A<-reshape(A,varying=list(c("age1","age2","age3","age4","age5","age6"),c("bmi1","bmi2","bmi3","bmi4","bmi5","bmi6"),
c("dbp1","dbp2","dbp3","dbp4","dbp5","dbp6")),idvar="id",
v.names=c("age","bmi","dbp"),times=c(1,2,3,4,5,6),direction="long")
names(A)
head(A)
tail(A)
A<-cc(A)   #complete data, na.omit(A) 
length(A$id)
# gender is a factor
A<-within(A, sex<-factor(Gender))


#############################################################
#Data analysis using package robustlmm
#############################################################
fit.robustlmer<-rlmer(dbp ~ age+sex+bmi+time+(1|id),data=A)