#Import data
data("PimaIndiansDiabetes2", package = "mlbench")
PimaIndiansDiabetes2

#Fig.1 (Missing values in the Pima Indian Diabetes dataset percentage bar plot and distribution.)
install.packages("VIM")
library(VIM)
mice_plot <- aggr(PimaIndiansDiabetes2, col=c('navyblue','yellow'),
                  numbers=TRUE, sortVars=TRUE,
                  labels=names(PimaIndiansDiabetes2), cex.axis=.7,
                  gap=3, ylab=c("Missing data","Pattern"))

#Fig.2 (Bivariate analysis with Pearson Correlation plot grouped by non-diabetic and diabetic patients. )
na.omit.data=na.omit(PimaIndiansDiabetes2)
library(ggplot2)
ggpairs(na.omit.data, columns = c("pregnant", "glucose", "pressure", "triceps","insulin","mass","pedigree","age"), title = "Bivariate analysis of revenue expenditure by the British household", upper = list(continuous = wrap("cor",
                                                                                                                                                                                                                                size = 3)),
        lower = list(continuous = wrap("smooth",
                                       alpha = 0.3,
                                       size = 0.1)),
        mapping = aes(color = diabetes))

#Fig.4 (The PCA eigenvalues and variances of the dimensions by bar and scree plot.)
#PCA visualize
library(SYNCSA)
pca.result.cluster5=pca(PimaIndiansDiabetes2[,1:8])
index.na.omit=as.numeric(rownames(na.omit(PimaIndiansDiabetes2)))
index.na=as.numeric(rownames(PimaIndiansDiabetes2)[apply(PimaIndiansDiabetes2, 1, anyNA)])
pca.result.cluster5$eigenvalues
# PCA plot1
library(ggplot2)
data <- data.frame(
  principal=c("1st","2nd","3rd","4th","5th","6th","7th","8th") ,  
  eigen.value=pca.result.cluster5$eigenvalues[,1]
)
ggplot(data=data, aes(x=principal, y=eigen.value)) +
  geom_bar(stat="identity", fill="steelblue")+
  geom_text(aes(label=round(eigen.value,3)), vjust=1.6, color="white", size=3.5)+
  theme_minimal()
# PCA plot2
data<-structure(list(principal=c("1st","2nd","3rd","4th","5th","6th","7th","8th"),percentage.intertia = pca.result.cluster5$eigenvalues[,2],cumulative.intertia = pca.result.cluster5$eigenvalues[,3]), .Names = c("principal", "percentage.intertia", "cumulative.intertia"), row.names = c(NA, 8L), class = "data.frame")
ggplot(data=data, aes(x=principal)) + 
  geom_bar(aes(y=percentage.intertia,fill="percentage intertia"), width=.7, stat="identity",fill="steelblue") +
  geom_line(aes(y=cumulative.intertia,group=1,linetype="Cumulative intertia"))+
  geom_point(aes(y=cumulative.intertia))+
  xlab("principal") + ylab("probability") +
  labs(fill="",linetype="")+
  theme_minimal()

#Fig. 5 (3-d scatter plot of three principal components by PCA with missing values for clustering to impute missing values.)
#PCA 3-dimension clustering plot
library(ClustImpute)
cluster.2=ClustImpute(PimaIndiansDiabetes2[,1:8],nr_cluster=4)
cluster.2$clusters
index.na.omit=as.numeric(rownames(na.omit(PimaIndiansDiabetes2)))
index.na=as.numeric(rownames(PimaIndiansDiabetes2)[apply(PimaIndiansDiabetes2, 1, anyNA)])
pch.vec=1:768
pch.vec[index.na.omit]=16
pch.vec[index.na]=17
library("scatterplot3d") # load
scatterplot3d(pca.result.cluster5$individuals[,1],pca.result.cluster5$individuals[,2],pca.result.cluster5$individuals[,3],pch=pch.vec,color = cluster.2$clusters+1,
              grid=TRUE, box=T ,scale.y = 2,xlab = "1st principal",ylab = "2nd principal",zlab = "3rd principal",angle = 10)
legend("bottom", legend = c("cluster 1", "cluster 1 (NA)","cluster 2","cluster 2 (NA)","cluster 3","cluster 3 (NA)","cluster 4","cluster 4 (NA)"),
       col =  c(2,2,3,3,4,4,5,5), 
       pch = c(16,17,16,17,16,17,16,17), 
       xpd = TRUE, horiz = TRUE, inset = -0.2,)

#Fig. 6
#Randomly create imputation test-validation dataset
install.packages("wakefield")
library(wakefield)
glucose.B=sample(x = c(T, F), 
                 prob = c(1-5/768, 5/768),
                 size = 392, 
                 replace = TRUE)
pressure.B=sample(x = c(T, F), 
                  prob = c(1-35/768, 35/768),
                  size = 392, 
                  replace = TRUE)
triceps.B=sample(x = c(T, F), 
                 prob = c(1-227/768, 227/768),
                 size = 392, 
                 replace = TRUE)
insulin.B=sample(x = c(T, F), 
                 prob = c(1-374/768, 374/768),
                 size = 392, 
                 replace = TRUE)
mass.B=sample(x = c(T, F), 
              prob = c(1-11/768, 11/768),
              size = 392, 
              replace = TRUE)
bool.df=!is.na(omit.PimaIndiansDiabetes.na)
bool.df[,2]=glucose.B
bool.df[,3]=pressure.B
bool.df[,4]=triceps.B
bool.df[,5]=insulin.B
bool.df[,6]=mass.B
omit.PimaIndiansDiabetes.na[!bool.df]=NA
#Middle: KNN imputation validation mean square error for K from 1 to 30 of weightAverage and median methods.
#distance measurement: "weighAvg"
knn.mse.weightA=c()
for(n.knn in 1:30){
  knn.n=KNNimp(omit.PimaIndiansDiabetes.na[,1:8], k = n.knn, scale = TRUE, meth = "weighAvg", distData = NULL)
  knn.mse.weightA=c(knn.mse.weightA,mean(((scale(knn.n)-scale(na.omit(PimaIndiansDiabetes2)[,-9]))[!bool.df[,-1]])^2))
  
}
knn.mse.weightA
#distance measurement: "median"
knn.mse.median=c()
for(n.knn in 1:30){
  knn.n=KNNimp(omit.PimaIndiansDiabetes.na[,1:8], k = n.knn, scale = TRUE, meth = "median", distData = NULL)
  knn.mse.median=c(knn.mse.median,mean(((scale(knn.n)-scale(na.omit(PimaIndiansDiabetes2)[,-9]))[!bool.df[,-1]])^2))
  
}
knn.mse.median
plot( 1:30,knn.mse.weightA , type="b" , bty="l" , xlab="k" , ylab="mean square error" , col=rgb(0.2,0.4,0.1,0.7) , lwd=3 , pch=17,
      main="KNN imputation validation-MSE for k from 1 to 30",cex.main=0.9)
lines(1:30,knn.mse.median , col=rgb(0.8,0.4,0.1,0.7) , lwd=3 , pch=19 , type="b" )
# Add a legend
legend("topright", 
       legend = c("weightAverage", "median"), 
       col = c(rgb(0.2,0.4,0.1,0.7), 
               rgb(0.8,0.4,0.1,0.7)), 
       pch = c(17,19), 
       bty = "n", 
       pt.cex = 2, 
       cex = 1.2, 
       text.col = "black", 
       horiz = F , 
       inset = c(0.1, 0.1))

#Right: PCA imputation validation mean square error for n-component from 1 to 7.
pca.mse=c()
for(n.pca in 1:7){
  pca.n=imputePCA(omit.PimaIndiansDiabetes.na[,-9],ncp=n.pca,scale = T)
  pca.mse=c(pca.mse,mean(((scale(pca.n$completeObs)-scale(na.omit(PimaIndiansDiabetes2)[,-9]))[!bool.df[,-1]])^2))
}
pca.mse
plot( 1:7,pca.mse , type="b" , bty="l" , xlab="number of components" , ylab="mean square error" , col=rgb(0.2,0.4,0.1,0.7) , lwd=3 , pch=17,
      main="PCA imputation validation-MSE for number of components from 1 to 7",cex.main=0.9)
#Left: Clustering imputation validation mean square error for clustering number of clusters from 1 to 15. 
dim(na.omit(PimaIndiansDiabetes2))
cluster.mse=c()
for(n.cluster in 1:15){
  cluster.n=ClustImpute(omit.PimaIndiansDiabetes.na[,1:8],nr_cluster=n.cluster)
  cluster.mse=c(cluster.mse,mean(((scale(cluster.n$complete_data)-scale(na.omit(PimaIndiansDiabetes2)[,-9]))[!bool.df[,-1]])^2))
  
}
plot( 1:15,cluster.mse , type="b" , bty="l" , xlab="number of clusters" , ylab="mean square error" , col=rgb(0.2,0.4,0.1,0.7) , lwd=3 , pch=17,
      main="clustering imputation validation-MSE for number of clusters from 1 to 15",cex.main=0.9)

#Fig. 7 (Multivariate normal distribution of the attributes with missing values and pedigree)
MCMC.impute.NA=PimaIndiansDiabetes2
MCMC.impute.NA$triceps=sqrt(MCMC.impute.NA$triceps)
MCMC.impute.NA$mass=sqrt(MCMC.impute.NA$mass)
MCMC.impute.NA$insulin=log(MCMC.impute.NA$insulin)
MCMC.impute.NA$pedigree=log(MCMC.impute.NA$pedigree)
names(MCMC.impute.NA)[4:7]=c("sqrt(triceps)","log(insulin)","sqrt(mass)","log(pedigree)")
library("PerformanceAnalytics")
chart.Correlation(MCMC.impute.NA[,2:7], histogram = T, pch= 19)

#MCMC impute missing values
mvn.data = MCMC.impute.NA[,2:7]
mvn.data
# Weakly informative prior parameter setting
n=dim(mvn.data)[1]
p<-dim(mvn.data)[2]
mu0<-apply(mvn.data,2,mean,na.rm=T)
sd0<-(mu0/2)
L0<-matrix(.1,p,p) ; diag(L0)<-1 ; L0<-L0*outer(sd0,sd0)
nu0<-p+2 ; S0<-L0
#Fully conditional posterior distribution for {θ,Σ,X_mis}:
#theta 
rmvnorm<-
  function(n,mu,Sigma) {
    p<-length(mu)
    res<-matrix(0,nrow=n,ncol=p)
    if( n>0 & p>0 ) {
      E<-matrix(rnorm(n*p),n,p)
      res<-t(  t(E%*%chol(Sigma)) +c(mu))
    }
    res
  }
theta.conditional.posterior = function(X.full, Sigma){
  ybar<-apply(X.full,2,mean)
  Ln<-solve( solve(L0) + n*solve(Sigma) )
  mun<-Ln%*%( solve(L0)%*%mu0 + n*solve(Sigma)%*%ybar )
  rmvnorm(1,mun,Ln)
}
#Sigma 
### sample from the Wishart distribution
rwish<-function(n,nu0,S0)
{
  sS0 <- chol(S0)
  S<-array( dim=c( dim(S0),n ) )
  for(i in 1:n)
  {
    Z <- matrix(rnorm(nu0 * dim(S0)[1]), nu0, dim(S0)[1]) %*% sS0
    S[,,i]<- t(Z)%*%Z
  }
  S[,,1:n]
}
Sigma.conditional.posterior = function(X.full,theta){
  Sn<- S0 + ( t(X.full)-c(theta) )%*%t( t(X.full)-c(theta) )
  solve( rwish(1, nu0+n, solve(Sn)) )
}
###missing matrix
O <-!is.na(mvn.data)
O<-O*1
THETA<-SIGMA<-X.MISS<-NULL
set.seed(1)
missing.index = which(apply(O,1,sum)<6)
#set starting values
Sigma<-S0
X.full<-mvn.data
for(j in 1:p)
{X.full[is.na(X.full[,j]),j]<-mean(X.full[,j],na.rm=TRUE)}
#MCMC iteration
S=5000
for(s in 1:S)
{
  ###update theta
  theta = theta.conditional.posterior(X.full,Sigma)
  ###update Sigma
  Sigma=Sigma.conditional.posterior(X.full,theta)
  ###update missing data
  for(i in missing.index){ 
    b <- ( O[i,]==0 )
    a <- ( O[i,]==1 )
    iSa<- solve(Sigma[a,a])
    beta.j <- Sigma[b,a]%*%iSa
    s2.j   <- Sigma[b,b] - Sigma[b,a]%*%iSa%*%Sigma[a,b]
    theta.j<- theta[b] + beta.j%*%(t(X.full[i,a])-theta[a])
    X.full[i,b] <- rmvnorm(1,theta.j,s2.j )
  }
  ### save results
  THETA<-rbind(THETA,theta) ; SIGMA<-rbind(SIGMA,c(Sigma))
  X.MISS<-rbind(X.MISS, X.full[O==0] )
  ###
  print(s)
}
thin = seq(1,5000,10)
#Fig.8(Convergence test plots (trace and acf) of sampling μ ̂_variable.)
par(mfrow=c(2,6))
plot(THETA[thin,1], main =expression(paste(theta[glucose])),xlab="iteration times/10",ylab = "value",type = "l",col="steelblue")
plot(THETA[thin,2],main=expression(paste(theta[pressure])),xlab="iteration times/10",ylab = "value",type = "l",col="steelblue")
plot(THETA[thin,3],main=expression(paste(theta[sqrt(triceps)])),xlab="iteration times/10",ylab = "value",type = "l",col="steelblue")
plot(THETA[thin,4],main=expression(paste(theta[log(insulin)])),xlab="iteration times/10",ylab = "value",type = "l",col="steelblue")
plot(THETA[thin,5],main=expression(paste(theta[sqrt(mass)])),xlab="iteration times/10",ylab = "value",type = "l",col="steelblue")
plot(THETA[thin,6],main=expression(paste(theta[log(pedigree)])),xlab="iteration times/10",ylab = "value",type = "l",col="steelblue")
acf(THETA[thin,1], main =expression(paste(theta[glucose])),xlab="Lag/10",col="steelblue")
acf(THETA[thin,2],main=expression(paste(theta[pressure])),xlab="Lag/10",col="steelblue")
acf(THETA[thin,3],main=expression(paste(theta[sqrt(triceps)])),xlab="Lag/10",col="steelblue")
acf(THETA[thin,4],main=expression(paste(theta[log(insulin)])),xlab="Lag/10",col="steelblue")
acf(THETA[thin,5],main=expression(paste(theta[sqrt(mass)])),xlab="Lag/10",col="steelblue")
acf(THETA[thin,6],main=expression(paste(theta[log(pedigree)])),xlab="Lag/10",col="steelblue")
X.missing.pred =apply(X.MISS[thin,],2,mean)
X.missing.pred

#Fig.9 (Visualization of the variability uncertainty on the plane defined by two PCA axes.)
install.packages("missMDA")
library(missMDA)
PimaIndiansDiabetes2
nb= estim_ncpPCA(PimaIndiansDiabetes2[,-9],scale = T)
comp=imputePCA(PimaIndiansDiabetes2[,-9],ncp=2,scale = T)
res.pca=PCA(comp$completeObs)
mi=MIPCA(PimaIndiansDiabetes2[,-9],scale = T,ncp=5)
plot(mi)

#mean and median method impute missing values result
install.packages("missMethods")
library(missMethods)
mean(((scale(impute_mean(omit.PimaIndiansDiabetes.na[,-9]))-scale(na.omit(PimaIndiansDiabetes2)[,-9]))[!bool.df[,-1]])^2)
mean(((scale(impute_median(omit.PimaIndiansDiabetes.na[,-9]))-scale(na.omit(PimaIndiansDiabetes2)[,-9]))[!bool.df[,-1]])^2)