# NMST 434: Exercise session 9

## Missing data, EM algorithm

rm(list=ls())
library(ellipse)
library(mvtnorm)

###
### 1. Missing completely at random 
###

# set.seed(1234567);
n <- 2000;
Sigma <- matrix(c(1,0.5,0.5,1), 2, 2);
X = rmvnorm(n,sigma=Sigma)

n0 <- 0.5*n;
n1 <- 0.25*n;
n2 <- 0.25*n;

# # Creating observed data 
Xo <- X;
Xo[(n0+1):(n0+n1),1] <- NA;
Xo[(n0+n1+1):(n0+n1+n2),2] <- NA;
Xo;

# in the following plot, each vertical line represents a (bivariate) vetor that represents a single observation. Its colour represents the value - darker color for higher values. Missing data are the white spaces.
image(Xo,xlab="Observations",ylab="Dimensions",main="Data value")

plotX = function(title = NULL){
  plot(X, cex = 0.5, xlab=expression(X[i]), ylab=expression(Y[i]), main = title); 
  points(0,0, col="blue", pch=4, cex=2, lwd=2)
  points(Xo, col="darkgreen", pch=8, cex=0.5)
  legend(2, -2, col=c("darkgreen", "black"),  pch=c(8,1), legend=c("observed", "missing"))
}

plotX()

# # A. Complete case analysis (CCA)
X.cc <- na.omit(Xo);  # complete cases 

(Sigma_CCA <- var(X.cc));
(mu_CCA <- apply(X.cc, 2, mean));

# Alternatively 
var(Xo);
var(Xo, na.rm = TRUE, use="complete.obs");

plotX("CCA")
points(mu_CCA[1],mu_CCA[2],col=2,pch=18,cex=2)
lines(ellipse(Sigma_CCA,centre = mu_CCA),col=2,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")

# # Confidence intervals for mu_1, mu_2
kvant <- qt(0.975, df=n0-1); 
CI_CCA <- matrix(mu_CCA, 2, 2) + kvant*matrix(c(-1,1), 2, 2, byrow=T)*sqrt(diag(Sigma_CCA))/sqrt(n0);
dimnames(CI_CCA) = list(c("mu_1","mu_2"),c("L","U"))
CI_CCA

# # B. Available case analysis (ACA)
# Variance matrix estimation 
(Sigma_ACA <- var(Xo, na.rm = TRUE, use="pairwise.complete"));

# To be compared with 
Sigma_CCA

# Alternatively 
var(Xo[,1], na.rm=T); 
var(Xo[,2], na.rm=T); 
cov(Xo[,1], Xo[,2], use="pairwise.complete"); 

# Mean estimation 
(mu_ACA <- apply(Xo, 2, mean, na.rm=TRUE));

# To be compared with 
mu_CCA

# Alternatively 
mean(na.omit(Xo[,1]))
mean(na.omit(Xo[,2]))

plotX("ACA")
points(mu_ACA[1],mu_ACA[2],col=2,pch=18,cex=2)
lines(ellipse(Sigma_ACA,centre = mu_ACA),col=2,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")

alfa <- 0.05;
# # Confidence interval for mu_1, mu_2
nr <- colSums(!is.na(Xo)); 
kvant <- qt(1-alfa/2, df=nr-1); 
CI_ACA <- matrix(mu_ACA, 2, 2) + matrix(c(-1,1), 2, 2, byrow=T) * matrix(kvant, 2, 2)*sqrt(diag(Sigma_CCA))/sqrt(nr);
dimnames(CI_ACA) = list(c("mu_1","mu_2"),c("L","U"))
CI_ACA

# # C. Imputation (simple)
# for imputation, we search for the best estimate of the missing data based on the observed data. For a bivariate normal model the conditional distribution of Y|X is normal again, so we can estimate E(Y|X) using a simple linear model
fit.12 <- lm(Xo[,1]~Xo[,2]);
X1hat <- predict(fit.12, data.frame(Xo[,2]));

fit.21 <- lm(Xo[,2]~Xo[,1]);
X2hat <- predict(fit.21, data.frame(Xo[,1])); 

# Filling data with imputed values 
X.imp <- Xo;
X.imp[is.na(X.imp[,1]),1]<- X1hat[is.na(X.imp[,1])]; 
X.imp[is.na(X.imp[,2]),2]<- X2hat[is.na(X.imp[,2])]; 

image(Xo,main="Observed data")
abline(v=c(.5,.75),col=3,lwd=4)
image(X.imp,main="Imputed data")
abline(v=c(.5,.75),col=3,lwd=4)
# you can see that the imputed data are not quite as the rest of the dataset..

(mu_imp <- apply(X.imp, 2, mean)); 
(Sigma_imp <- var(X.imp));

# To be compared with 
mu_ACA
Sigma_ACA

plotX("Simple Imputation")
points(mu_imp[1],mu_imp[2],col=2,pch=18,cex=2)
lines(ellipse(Sigma_imp,centre = mu_imp),col=2,lwd=2)
lines(ellipse(Sigma_ACA,centre = mu_ACA),col=3,lwd=2)
lines(ellipse(Sigma_CCA,centre = mu_CCA),col=4,lwd=2)
lines(ellipse(Sigma),col=1,lwd=3, lty="dotted")
legend("topleft",c("Simple Imputation","ACA","CCA","True"),col=c(2,3,4,1),lwd=c(2,2,2,3), lty=c(rep("solid", 3),"dotted"))

# plot of the imputed values, vertical lines represent where only the X-coordinate of a point is known. At this line, the red point is the  imputed value. The orange line is the estimated regression line Y~X
plotX("Simple Imputation")
abline(v=Xo[is.na(Xo[,2]),1],lwd=.05,lty=3)
abline(fit.21,col="goldenrod",lwd=3)
points(Xo[is.na(Xo[,2]),1],X2hat[is.na(Xo[,2])],pch=16,col=2)

# jointly, all imputed values lie only on the two estimated regression 
# lines Y~X and X~Y
plotX("Simple Imputation")
points(X.imp[is.na(Xo[,2]),1],X.imp[is.na(Xo[,2]),2],pch=16,col=2,cex=.5)
points(X.imp[is.na(Xo[,1]),1],X.imp[is.na(Xo[,1]),2],pch=16,col=2,cex=.5)

# # Confidence interval for mu_1, mu_2 - what is the problem here?
(CI_imp <- matrix(mu_imp, 2, 2) + matrix(c(-1,1), 2, 2, byrow=T) * qt(0.975, df=n-1) * sqrt(diag(Sigma_imp))/sqrt(n));


# # D. Multiple imputation 
# We now assume that the data is indeed bivariate normal, and estimate not only the conditional expectation of the missing data given the observed data, but also its variability and consider this in our computation
sigma.12 <- sqrt(sum(fit.12$residuals^2)/(n0-2));
sigma.21 <- sqrt(sum(fit.21$residuals^2)/(n0-2));

# Filling data with imputed values 
na.X <- sum(is.na(Xo[,1]));
na.Y <- sum(is.na(Xo[,2]));

M <- 1000; # Number of imputations 
mu.imp.M <- matrix(NA, nrow=M, ncol=2);
Sigma.imp.M <- matrix(NA, nrow=M, ncol=4); 

# we repeat the imputation procedure independently M times. 
# in each run, we estimate mu and Sigma. This will give us 
# some information about the variability of the estimator due
# to imputing the missing data
for(i in 1:M){
  X.imp <- Xo;
  X.imp[is.na(X.imp[,1]),1] <- X1hat[is.na(X.imp[,1])] + sigma.12 * rnorm(na.X); 
  X.imp[is.na(X.imp[,2]),2] <- X2hat[is.na(X.imp[,2])] + sigma.21 * rnorm(na.Y); 

  mu.imp.M[i,] <- apply(X.imp, 2, mean); 
  Sigma.imp.M[i,] <- as.vector(var(X.imp));
}


plotX("Multiple Imputation")
points(X.imp[is.na(Xo[,2]),1],X.imp[is.na(Xo[,2]),2],pch=3,col="blue",cex=0.75)
points(mu.imp.M[M,1],mu.imp.M[M,2],col=2,pch=18,cex=2)
lines(ellipse(matrix(Sigma.imp.M[M,],nrow=2),centre = mu.imp.M[M,]),col=2,lwd=2)
lines(ellipse(Sigma), col=1,lwd=3, lty="dotted")
legend("topleft",c("Parametric Imputation","True"),col=c(2,1),lwd=c(2,3), lty=c("solid", "dotted"))

# aggregated estimator from all M imputations
(mu_imp_mult <- colMeans(mu.imp.M));
(Sigma_imp_mult <- matrix(colMeans(Sigma.imp.M), 2, 2));

# To be compared with 
mu_imp
Sigma_imp;


# Q: What can be said about mu_imp_mult in comparison with mu_imp ? 

# # Confidence interval construction (see Section 7.2 in the course notes)
# Variance estimation of mu_imp_mult 
(VM <- colMeans(Sigma.imp.M[,c(1,4)])/n); 
# Estimation of expectation of conditional variance 
(BM <- diag(var(mu.imp.M)));  
# Estimation of variance of conditional expectation 
(var.mult.imp <- VM + BM);
(CI_imp_mult <- matrix(mu_imp, 2, 2) + matrix(c(-1,1), 2, 2, byrow=T) * 
   matrix(qt(0.975, df=nr-1), 2, 2) * sqrt(var.mult.imp));

# Think about the degrees of freedom used in the quantile of t-distribution. 

# To be compared with
CI_imp;


# # E. EM - algorithm 

# Initial estimates 
mu_EM <- mu_CCA; 
Sigma_EM <- Sigma_CCA; 

n.iter <- 100;

# Observed log-likelihood 
loglik <- function(n0, n1, n2, mu, Sigma){
  Mu_hat <- matrix(mu, n0, 2, byrow=TRUE); 
  -n0*log(det(Sigma)) - n1*log(Sigma[2,2]) - n2*log(Sigma[1,1]) - 
    sum(diag((Xo[1:n0,]-Mu_hat)%*%solve(Sigma)%*%t((Xo[1:n0,]-Mu_hat)))) - 
    sum((Xo[(n0+n1+1):(n0+n1+n2),1] - mu[1])^2)/Sigma[1,1] - 
    sum((Xo[(n0+1):(n0+n1),2] - mu[2])^2)/Sigma[2,2];
}

loglik.old <-loglik(n0, n1, n2, mu_EM, Sigma_EM); 

# loglik.old
# library(mvtnorm)
# sum(dmvnorm(na.omit(Xo),mean=mu_EM, sigma=Sigma_EM,log=TRUE))+
#   sum(dnorm(Xo[is.na(Xo[,2]),1],mean=mu_EM[1],sd=sqrt(Sigma_EM[1,1]),log=TRUE))+
#   sum(dnorm(Xo[is.na(Xo[,1]),2],mean=mu_EM[2],sd=sqrt(Sigma_EM[2,2]),log=TRUE))
# (loglik.old-(2*n0+n1+n2)*log(2*pi))/2

for(i in 1:n.iter){
  
  beta <- Sigma_EM[1,2]/c(Sigma_EM[2,2], Sigma_EM[1,1]);
  Yhat <- mu_EM[2] + beta[2]*(Xo[,1] - mu_EM[1]);  # predicted Y
  Xhat <- mu_EM[1] + beta[1]*(Xo[,2] - mu_EM[2]);  # predicted X 

  rho <- (Sigma_EM[1,2])/sqrt(Sigma_EM[1,1]*Sigma_EM[2,2]); # correlation coefficient 
    
  Yhat2 <- Sigma_EM[2,2]*(1 - rho^2) + Yhat^2; # predicted Y^2 
  Xhat2 <- Sigma_EM[1,1]*(1 - rho^2) + Xhat^2; # predicted X^2
  
  XX <- Xo;
  XX[(n0+1):(n0+n1),1] <- Xhat[(n0+1):(n0+n1)];
  XX[(n0+n1+1):(n0+n1+n2),2] <- Yhat[(n0+n1+1):(n0+n1+n2)];

  mu_EM <- colMeans(XX);
  
  XX2 <- Xo^2;
  XX2[(n0+1):(n0+n1),1] <- Xhat2[(n0+1):(n0+n1)];
  XX2[(n0+n1+1):(n0+n1+n2),2] <- Yhat2[(n0+n1+1):(n0+n1+n2)];

  diag(Sigma_EM) <- colMeans(XX2) - mu_EM^2;
  Sigma_EM[1,2] <- mean(XX[,1]*XX[,2]) - mu_EM[1] * mu_EM[2];
  Sigma_EM[2,1] <- Sigma_EM[1,2];
  
  loglik.new <- loglik(n0, n1, n2, mu_EM, Sigma_EM); 
  
  print(c(i, loglik.new));

# #   Stopping rule 
  if(loglik.new-loglik.old < 1e-10){
    break;
  }else{
    loglik.old <- loglik.new;
  }  
}

print(mu_EM);
print(Sigma_EM);

# With the help of library 'norm2'
# install.packages(norm2)

library(norm2);
summary(emNorm(Xo))


# the same using direct optimization
loglik <- function(th){
  Mu_hat <- matrix(th[1:2], n0, 2, byrow=TRUE);
  Sigma = matrix(c(th[3],th[4],th[4],th[5]),ncol=2)
  -(-n0*log(det(Sigma)) - n1*log(Sigma[2,2]) - n2*log(Sigma[1,1]) - 
    sum(diag((Xo[1:n0,]-Mu_hat)%*%solve(Sigma)%*%t((Xo[1:n0,]-Mu_hat)))) - 
    sum((Xo[(n0+n1+1):(n0+n1+n2),1] - th[1])^2)/Sigma[1,1] - 
    sum((Xo[(n0+1):(n0+n1),2] - th[2])^2)/Sigma[2,2]);
}

opt = optim(c(0,0,1,1/2,1),loglik)
opt$par[1:2]
matrix(opt$par[c(3,4,4,5)],ncol=2)

# comparison of objective functions (should be maximized)
-opt$value
loglik.new
2*summary(emNorm(Xo))$loglik

# # # Summary 

# Estimation of parameters 
mu_res <- rbind(c(mu_CCA, as.numeric(Sigma_CCA)[-2]), c(mu_ACA, as.numeric(Sigma_ACA)[-2]), 
	c(mu_imp, as.numeric(Sigma_imp)[-2]), c(mu_imp_mult, as.numeric(Sigma_imp_mult)[-2]), 
	c(mu_EM, as.numeric(Sigma_EM)[-2]));
colnames(mu_res) <- c("E(X)", "E(Y)", "var(X)", "cov(X,Y)", "var(Y)");
rownames(mu_res) <- c("CCA", "ACA", "imp", "imp_mult", "mu_EM");
print(mu_res, 2) 


# Confidence interval for mu_{2}
CI_mu2 <- rbind(CI_CCA[2,], CI_ACA[2,], CI_imp[2,], CI_imp_mult[2,]);
CI_mu2 = cbind(CI_mu2,CI_mu2[,2]-CI_mu2[,1])
colnames(CI_mu2) <- c("lower", "upper","length");
rownames(CI_mu2) <- c("CCA", "ACA", "imp", "imp_mult");
print(CI_mu2, 2) 

# Q: How to construct the CI with the help of "ignorable" likelihood method? 


### 
### 2. Missing at random   
### 

source("fctions_missing.R") 
# some elementary functions from the above, stored in an external file for 
# convenience

# set.seed(1234567)
n <- 1000;
X = rmvnorm(n,sigma=Sigma)

# # Creating observed data 
# censoring function: if the values is too small, the datum will not be observed
cen.fun = function(x) exp(4*x)/(1+exp(4*x))
# cen.fun = function(x) dnorm(x,sd=.5)*sqrt(2*pi)*.5

Xo <- X;
qqq <- as.logical(cen.fun(X[,1]) < runif(n));
Xo[qqq,2] <- NA;
Xo
summary(Xo)
image(Xo)

plotX()

# 1. Complete case analysis 
(res.CCA <- binormal.CCA(Xo))

plotX("CCA")
points(res.CCA$muhat[1],res.CCA$muhat[2],col=2,pch=18,cex=2)
lines(ellipse(res.CCA$Sigma2hat,centre = res.CCA$muhat),col=2,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")


# 2. Available case analysis 
(res.ACA <- binormal.ACA(Xo))

plotX("ACA")
points(res.ACA$muhat[1],res.ACA$muhat[2],col=2,pch=18,cex=2)
lines(ellipse(res.ACA$Sigma2hat,centre = res.ACA$muhat),col=2,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")

# 3. Imputation
(res.imp <- binormal.imp(Xo))

plotX("Simple Imputation")
points(res.imp$muhat[1],res.imp$muhat[2],col=2,pch=18,cex=2)
lines(ellipse(res.imp$Sigma2hat,centre = res.imp$muhat),col=2,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")

# 4. Multiple imputation
(res.mult.imp <- binormal.imp.mult(Xo))

plotX("Multiple Imputation")
points(res.mult.imp$muhat[1],res.mult.imp$muhat[2],col=2,pch=18,cex=2)
lines(ellipse(res.mult.imp$Sigma2hat,centre = res.mult.imp$muhat),col=2,lwd=2)
lines(ellipse(Sigma_imp,centre = mu_imp),col=3,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")
legend("topleft",c("Multiple Imputation","Simple Imputation","True"),col=c(2,3,1),lwd=c(2,2,3), lty=c("solid", "solid", "dotted"))

# 5. EM algorithm
(res.EM <- binormal.EM(Xo))

plotX("EM algorithm")
points(res.EM$muhat[1],res.EM$muhat[2],col=2,pch=18,cex=2)
lines(ellipse(res.EM$Sigma2hat,centre = res.EM$muhat),col=2,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")

# 6. Reweighting 

# Estimating EY and var(Y)
R <- as.logical(!is.na(Xo[,2])); # indicator of the observed values 
(fit <- glm(R~Xo[,1], family=binomial));
summary(pihat <- fit$fitted);
plot(Xo[,1], pihat, xlab=expression(X[i]), ylab=expression(pi[i]),ylim=c(0,1))
points(Xo[,1],cen.fun(Xo[,1]),cex=.5,col=2)
legend("topleft",c("Estimated w[i]", "True w[i]"),col=1:2,pch=1)

wi <- (1/pihat);
wi[qqq] <- NA;
wi <- wi/sum(wi, na.rm=T)
(mu2_rew <- sum(Xo[,2]*wi, na.rm=T))
(sigma2_rew <- sum((Xo[,2]-mu2_rew)^2*wi, na.rm=T))

# plot(Xo[!is.na(wi),1], wi[!is.na(wi)], xlab=expression(X[i]), ylab=expression(w[i]))

(sigma12_rew <- sum((Xo[,1]-res.ACA$muhat[,1])*(Xo[,2]-mu2_rew)*wi, na.rm=T)/(sqrt(sigma2_rew * res.ACA$Sigma2hat[1,1])))

# Estimated Sigma matrix
Xoc = t(t(Xo)-c(res.ACA$muhat[1],mu2_rew))[!qqq,] # complete observations centred at the estimated mean vector
(Sigma2hat = t(wi[!qqq]*Xoc)%*%Xoc)

# cov2cor(Sigma2hat) # not the same off-diagonal term as sigma12_rew. Why?

plotX("Reweighting")
points(res.ACA$muhat[1],mu2_rew,col=2,pch=18,cex=2)
lines(ellipse(Sigma2hat,centre = c(res.ACA$muhat[1],mu2_rew)),col=2,lwd=2)
# 
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")
legend("topleft",c("Reweighting","True"),col=c(2,1), lwd=c(2,3), lty=c("solid", "dotted"))

# # Summary 
# Estimation of parameters 
mu_res <- rbind(c(res.CCA$muhat[,1], as.numeric(res.CCA$Sigma2hat)[-2]), c(res.ACA$muhat[,1], as.numeric(res.ACA$Sigma2hat)[-2]), c(res.imp$muhat, as.numeric(res.imp$Sigma2hat)[-2]), 
	c(res.mult.imp$muhat[,1], as.numeric(res.mult.imp$Sigma2hat)[-2]), c(res.EM$muhat, as.numeric(res.EM$Sigma2hat)[-2]), c(mean(X[,1]), mu2_rew, var(X[,1]), sigma12_rew, sigma2_rew));
colnames(mu_res) <- c("E(X)", "E(Y)", "var(X)", "cov(X,Y)", "var(Y)");
rownames(mu_res) <- c("CCA", "ACA", "imp", "imp_mult", "mu_EM", "reweight.");
print(mu_res, 2) 

# Confidence interval for mu_{1}
CI_mu1 <- rbind(res.CCA$muhat[1,2:3], res.ACA$muhat[1,2:3], res.mult.imp$muhat[1,2:3]);
colnames(CI_mu1) <- c("lower", "upper");
rownames(CI_mu1) <- c("CCA", "ACA", "imp_mult");
print(CI_mu1, 2) 


# # Confidence interval for mu_{2}

CI_mu2 <- rbind(res.CCA$muhat[2,2:3], res.ACA$muhat[2,2:3], res.mult.imp$muhat[2,2:3]);
CI_mu2 = cbind(CI_mu2, CI_mu2[,2]-CI_mu2[,1])
colnames(CI_mu2) <- c("lower", "upper","length");
rownames(CI_mu2) <- c("CCA", "ACA", "imp_mult");
print(CI_mu2, 2) 

# # # # # # # # # # # # # 
# # Test the null hypothesis that mu0 = mu1
X.cc <- na.omit(Xo); # complete cases 

t.test(X.cc[,1]-X.cc[,2]) # t-test for complete cases only

(loglik1 <- binormal.EM(Xo)$loglik); 
summary(emNorm(Xo))

(loglik0 <- EM2(Xo)$loglik); # log-likelihood under the null hypothesis 

# p-value of the likelihood ratio test
1 - pchisq(2*(loglik1-loglik0), df=1)

### 
### 3. Missing not at random 
### 

n <- 1000;
X = rmvnorm(n,sigma=Sigma)

# # Creating observed data 

cen.fun = function(x,y) exp(x-y)/(1+exp(x-y))

Xo <- X;
qqq <- as.logical(cen.fun(X[,1],X[,2]) < 0.5);
Xo[qqq,2] <- NA;
Xo;
summary(Xo)
image(Xo)

plotX()
contour(x<-sort(X[,1]),y<-sort(X[,2]),outer(x,y,cen.fun),add=TRUE,levels=1/2)

# 1. Complete case analysis 
(res.CCA <- binormal.CCA(Xo))

plotX("CCA")
points(res.CCA$muhat[1],res.CCA$muhat[2],col=2,pch=18,cex=2)
lines(ellipse(res.CCA$Sigma2hat,centre = res.CCA$muhat),col=2,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")

# 2. Available case analysis 
(res.ACA <- binormal.ACA(Xo))

plotX("ACA")
points(res.ACA$muhat[1],res.ACA$muhat[2],col=2,pch=18,cex=2)
lines(ellipse(res.ACA$Sigma2hat,centre = res.ACA$muhat),col=2,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")

# 3. Imputation
(res.imp <- binormal.imp(Xo))

plotX("Imputation")
points(res.imp$muhat[1],res.imp$muhat[2],col=2,pch=18,cex=2)
lines(ellipse(res.imp$Sigma2hat,centre = res.imp$muhat),col=2,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")

# 4. Multiple imputation
(res.mult.imp <- binormal.imp.mult(Xo))

plotX("Multiple Imputation")
points(res.mult.imp$muhat[1],res.mult.imp$muhat[2],col=2,pch=18,cex=2)
lines(ellipse(res.mult.imp$Sigma2hat,centre = res.mult.imp$muhat),col=2,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")

# 5. EM algorithm
(res.EM <- binormal.EM(Xo))

plotX("EM Algorithm")
points(res.EM$muhat[1],res.EM$muhat[2],col=2,pch=18,cex=2)
lines(ellipse(res.EM$Sigma2hat,centre = res.EM$muhat),col=2,lwd=2)
lines(ellipse(Sigma,centre = mu_CCA),lwd=3, lty="dotted")

# 6. Reweighting 
R <- as.logical(!is.na(Xo[,2])); # indicatior of the observed values 
(fit <- glm(R~Xo[,1], family=binomial));
summary(pihat <- fit$fitted);

plot(Xo[,1], pihat, xlab=expression(X[i]), ylab=expression(pi[i]),ylim=c(0,1))
points(X[,1],cen.fun(X[,1],X[,2]),cex=.5,col=2)
abline(h=.5,col=2)
legend("topleft",c("Estimated w[i]", "Censoring function value","Cutoff for censoring"),col=c(1,2,2),pch=c(1,1,NA),lty=c(0,0,1))

wi <- (1/pihat);
wi[qqq] <- NA;
wi <- wi/sum(wi, na.rm=T)
(mu2_rew <- sum(Xo[,2]*wi, na.rm=T))
(sigma2_rew <- sum((Xo[,2]-mu2_rew)^2*wi, na.rm=T))

# Estimated Sigma matrix
Xoc = t(t(Xo)-c(res.ACA$muhat[1],mu2_rew))[!qqq,] # complete observations centred at the estimated mean vector
(Sigma2hat = t(wi[!qqq]*Xoc)%*%Xoc)

plotX("Reweighting")
points(res.ACA$muhat[1],mu2_rew,col=2,pch=18,cex=2)
lines(ellipse(Sigma2hat,centre = c(res.ACA$muhat[1],mu2_rew)),col=2,lwd=2)
#
lines(ellipse(Sigma), col=1,lwd=3, lty="dotted")
legend("topleft",c("Reweighting","True"),col=c(2,1),lwd=c(2,3), lty=c("solid", "dotted"))


# # Summary 
# Estimation of parameters 
mu_res <- rbind(c(res.CCA$muhat[,1], as.numeric(res.CCA$Sigma2hat)[-2]), c(res.ACA$muhat[,1], as.numeric(res.ACA$Sigma2hat)[-2]), c(res.imp$muhat, as.numeric(res.imp$Sigma2hat)[-2]), 
	c(res.mult.imp$muhat[,1], as.numeric(res.mult.imp$Sigma2hat)[-2]), c(res.EM$muhat, as.numeric(res.EM$Sigma2hat)[-2]), c(NA, mu2_rew, NA, NA, sigma2_rew));
colnames(mu_res) <- c("E(X)", "E(Y)", "var(X)", "cov(X,Y)", "var(Y)");
rownames(mu_res) <- c("CCA", "ACA", "imp", "imp_mult", "mu_EM", "reweight.");
print(mu_res, 2) 

# Confidence interval for mu_{1}
CI_mu1 <- rbind(res.CCA$muhat[1,2:3], res.ACA$muhat[1,2:3], res.mult.imp$muhat[1,2:3]);
CI_mu1 = cbind(CI_mu1,CI_mu1[,2]-CI_mu1[,1])
colnames(CI_mu1) <- c("lower", "upper","length");
rownames(CI_mu1) <- c("CCA", "ACA", "imp_mult");
print(CI_mu1, 2) 

# Confidence interval for mu_{2}
CI_mu2 <- rbind(res.CCA$muhat[2,2:3], res.ACA$muhat[2,2:3], res.mult.imp$muhat[2,2:3]);
CI_mu2 = cbind(CI_mu2,CI_mu2[,2]-CI_mu2[,1])
colnames(CI_mu2) <- c("lower", "upper","length");
rownames(CI_mu2) <- c("CCA", "ACA", "imp_mult");
print(CI_mu2, 2) 

###
### 4. Linear regression with missing values in a covariate 
### 

n <- 1000;
X = rmvnorm(n,sigma=Sigma) # normally distributed covariates 
# X <- 2*sqrt(3)*(matrix(runif(2*n), 2, n) - 0.5); # uniformly distributed covariates

beta <- c(0,1)
Y <- X%*%beta + rnorm(n);

# # Creating observed data 
Xo <- X;
qqq <- as.logical(exp(X[,1]+Y)/(1+exp(X[,1]+Y)) < 0.5);
Xo[qqq,2] <- NA;
summary(Xo)
image(cbind(Xo,Y))

plot(X, cex = 0.5, xlab=expression(X[1]), ylab=expression(X[2])); 
points(Xo, col="darkgreen", pch=8, cex=0.5)
legend("bottomright", col=c("darkgreen", "black"),  pch=c(8,1), legend=c("observed", "missing"))

apply(cbind(Xo,Y), 2, mean, na.rm=T)

summary(fit.true <- lm(Y~X))

# CCA 
summary(fit.CCA <- lm(Y~Xo));
(CI.CCA <- confint(fit.CCA))

# Direct likelihood approach assuming trivariate normality of data 
temp <- summary(emNorm(cbind(Xo,Y)));
V <- temp$param$sigma;
betahat_slope_EM <- solve(V[1:2,1:2], V[1:2,3]);
(betahat_EM <- c(temp$param$beta[3] - betahat_slope_EM%*%temp$param$beta[1:2], betahat_slope_EM))

# Multiple imputation method 
fit <- lm(Xo[,2]~Y + Xo[,1]);
sigmahat <- sqrt(sum(fit$res^2)/fit$df); 

X2hat <- cbind(1, Y, Xo[,1])%*%fit$coeff; 

M <- 20; 
X.imp <- Xo;
qqq <- is.na(Xo[,2]);
n.m <- sum(qqq)
beta.M <- matrix(NA, nrow=M, ncol=3); 
var.M <- matrix(NA, nrow=M, ncol=3); 

for(i in 1:M){  
  X.imp[qqq,2] <- X2hat[qqq] + sigmahat*rnorm(n.m); 
  fitM <- lm(Y~X.imp);
  beta.M[i,] <- fitM$coef;
  var.M[i,] <- (summary(fitM)$coeff[,2])^2;
}

# summary(beta.M);
# summary(var.M);

# Estimate of parameters beta 
(beta.MI <- colMeans(beta.M));
(var.MI <- colMeans(var.M) + (1+1/M)*apply(beta.M, 2, var));

beta_res <- rbind(fit.CCA$coef, beta.MI, betahat_EM, fit.true$coef);
colnames(beta_res) <- c("beta0", "beta1", "beta2");
rownames(beta_res) <- c("CCA", "imp_beta", "beta_EM", "Complete");
print(beta_res, 2) 

# # Confidence intervals 

# Multiple imputation 
CI.imp <- matrix(beta.MI, nrow=3, ncol=2) + matrix(c(-1,1), nrow=3, ncol=2, byrow=T) * qt(0.975, df=n-n.m-3)*matrix(sqrt(var.MI), nrow=3, ncol=2);
colnames(CI.imp) <- c("lower", "upper");
rownames(CI.imp) <- c("intercept", "X1", "X2");
print(CI.imp, digits=2)


# Beta_0
CI.beta0 <- rbind(CI.CCA[1,], CI.imp[1,], confint(fit.true)[1,]);
rownames(CI.beta0) <- c("CCA", "Mult. imp", "Complete")
print("Beta_0:")
print(CI.beta0, digits=2)

# # Beta_0
# CI.beta0 <- rbind(CI.CCA[1,], CI.imp[1,], confint(fit.true)[1,]);
# rownames(CI.beta0) <- c("CCA", "Mult. imp", "Complete")
# print("Beta_0:")
# print(CI.beta0, digits=2)

# Beta_1
CI.beta1 <- rbind(CI.CCA[2,], CI.imp[2,], confint(fit.true)[2,]);
rownames(CI.beta1) <- c("CCA", "Mult. imp", "Complete")
print("Beta_1:")
print(CI.beta1, digits=2)

# Beta_2
CI.beta2 <- rbind(CI.CCA[3,], CI.imp[3,], confint(fit.true)[3,]);
rownames(CI.beta2) <- c("CCA", "Mult. imp", "Complete")
print("Beta_2:")
print(CI.beta2, digits=2)