###
###  MODERN STATISTICAL METHODS (NMST434)
###  Exercise class -- week 5 
###

rm(list=ls());

# # # # # # # # # # # # # # # # # # # # # # # # # # # #  
# # # 1. Misspecified poisson regression 
# # # # # # # # # # # # # # # # # # # # # # # # # # # # 

load("data/concentration_vs_counts.RData")

# The dataset DATA contains the results of 100 observations.  Each of the observation  contains information about a concentration of a given protein in the sample and the the count (number) of bacteries found in the sample. The task is to describe how the count is influenced by the concentration. 

summary(DATA)
plot(DATA$concentration, DATA$count, xlab="concentration", ylab="count")

summary(fit <- glm(count~concentration, data=DATA, family="poisson"))

X <- model.matrix(fit);
n <- nrow(X);
Y <- DATA$count;

# The sandwich estimate of the asymptotic variance given at the lecture can be calculated as follows. 
# Estimate of the matrix Gamma 
(Gammahat <- (1/n)*t(X)%*%diag(fitted(fit))%*%X);  
# Estimate of the matrix Sigma  
(Sigmahat <- (1/n)*t(X)%*%diag(as.vector(Y-fitted(fit))^2)%*%X);   

# When believing fully our Poisson regression model 
(1/n)*solve(Gammahat); 

# Compare with the estimate of asymptotic variance as given by glm(, family="poisson") function 
vcov(fit)
vcov(fit) - (1/n)*solve(Gammahat)

# Sandwich estimator of the asymptotic variance 
(estim.avar <- (1/n)*solve(Gammahat)%*%Sigmahat%*%solve(Gammahat)); 

# Note that this estimator coincides with the sandwich estimator from the library 'sandwich'
library(sandwich);
sandwich(fit);
sandwich(fit) - estim.avar;


                                         
# # # Small simulation experiment 
nopak <- 100;
beta <- c(1,1);
CI <- matrix(NA, nopak, 2)
CI.sand <- matrix(NA, nopak, 2)
alfa <- 0.05;
kvant <- qnorm(1-alfa/2);
beta <- c(1,1);

for(i in 1:nopak){
	print(i)
	n <- 100;  # sample size 
	x0 <- rep(1, n);        # regressor for the intercept  
	x1 <- floor(9*runif(n))/2;       # regressor for the slope 
	# regressor for the slope
	X <- cbind(x0, x1);     # the regression matrix X
	
	# 	Model misspecification, The conditional distribituion of Y|X is Bin(2, 0.5)*exp(x*beta)
	Y <- exp(X%*%beta)*(rbinom(n, size=20, prob=0.5)/10);
	
	fit <- glm(Y~x1, family="poisson");
	beta1 <- coef(fit)[2];
	CI[i,] <- beta1 + c(-1,1) * kvant*sqrt(vcov(fit)[2,2]);
	CI.sand[i,] <- beta1 + c(-1,1) * kvant*sqrt(sandwich(fit)[2,2]);	
# 	plot(x1, Y)
}

# # Assessing empirical properties of CI
assess.CI <- function(CI, theta){
	coverage <- mean((CI[,1] <= theta)*(theta <= CI[,2]));
	aver.length <- mean(CI[,2]-CI[,1]);
	sd.length <- sd(CI[,2]-CI[,1]);
  
	ret <- c(coverage, aver.length, sd.length);
	names(ret) <- c("coverage", "mean(length)", "sd(length)");
	
	ret;
}


# Confidence intervals for beta_1 (i.e. beta[2])
res.poiss <- assess.CI(CI, beta[2]);
res.sand <- assess.CI(CI.sand, beta[2]);

round(rbind(res.poiss, res.sand),3)

# # # # Model misspecification, n=100, nopak=10000
#           coverage mean(length) sd(length)
# res.poiss    0.783        0.071      0.003
# res.sand     0.926        0.108      0.020



# # # # Model misspecification, n=500, nopak=10000
#           coverage mean(length) sd(length)
# res.poiss    0.776        0.032      0.001
# res.sand     0.946        0.050      0.004




# # # # # # # # # # # # # # # # # # # # # # # # # # # #                                          
# 2. White estimator of asymptotic variance in a linear model 
# # # # # # # # # # # # # # # # # # # # # # # # # # # #                                          

library(sandwich);

n <- 100;
x0 <- rep(1, n);        # regressor for the intercept  
x1 <- 2*runif(n);       # regressor for the slope
X <- cbind(x0, x1);     # the regression matrix X

beta <- c(1,1);         # true values of the parameters  

# # Heteroscedastic normal errors 
sigma <- function(x1) exp(x1)
eps <- (rexp(n)-1)*sigma(x1);
Y <- X%*%beta + eps;

# Summary 
summary(fit.lm <- lm(Y~x1));

# Plot 
plot(x1, Y, cex=0.2)
abline(fit.lm$coef, col="blue", lwd=2)

# Sandwich estimator (White's estimator)
sandwich(fit.lm); 

# What is being computed ? 
# Estimate of the matrix Gamma 
(Gammahat <- (1/n)*t(X)%*%X);
# Estimate of the matrix Sigma  
(Sigmahat <- (1/n)*t(residuals(fit.lm)^2*X)%*%X);   

# Sandwich estimator of the asymptotic variance - textbook version  
(estim.avar <- (1/n)*solve(Gammahat)%*%Sigmahat%*%solve(Gammahat)); 

# To improve finite sample performance of the variance 
# it is recommended that the squared residuals are 
# divided by (1 - h_{ii})^2
vcovHC(fit.lm, type = "HC3")

# Manual calculation 
# Estimate of the Sigma matrix - version HC3
(Sigmahat.HC3 <- (1/n)*t(residuals(fit.lm)^2/(1-hatvalues(fit.lm))^2*X)%*%X);   
# Sandwich estimate of the variance matrix - version HC3
(1/n)*solve(Gammahat)%*%Sigmahat.HC3%*%solve(Gammahat);

# vcovHC calculation minus manual calculation 
vcovHC(fit.lm, type = "HC3") - (1/n)*solve(Gammahat)%*%Sigmahat.HC3%*%solve(Gammahat);

# Comparison 
summary(fit.lm);
sqrt(diag(sandwich(fit.lm)))
sqrt(diag(vcovHC(fit.lm, type = "HC3")))

# # Small simluation study 
# Comparison of the different estimators of asymptotic variance through 
# confidence intervals 

nopak <- 100;
n <- 100;
beta1 <- 1;

# Confidence intervals intercept 
CI.stand <- matrix(NA, nrow=nopak, ncol=2);
CI.white <- matrix(NA, nrow=nopak, ncol=2);
CI.HC3 <- matrix(NA, nrow=nopak, ncol=2);

# Confidence intervals slopes 
CIs.stand <- matrix(NA, nrow=nopak, ncol=2);
CIs.white <- matrix(NA, nrow=nopak, ncol=2);
CIs.HC3 <- matrix(NA, nrow=nopak, ncol=2);

# Variance function 
sigma <- function(x1) exp(x1)

kvant <- qt(0.975, df=n-2);
# # # Assessing empirical properties of CI
# assess.CI <- function(CI, theta){
# 	coverage <- mean((CI[,1] <= theta)*(theta <= CI[,2]));
# 	aver.length <- mean(CI[,2]-CI[,1]);
# 	sd.length <- sd(CI[,2]-CI[,1]);
#   
# 	ret <- c(coverage, aver.length, sd.length);
# 	names(ret) <- c("coverage", "mean(length)", "sd(length)");
# 	
# 	ret;
# }

for(i in 1:nopak){
  print(i);
# # 	Data generation process 
	x1 <- 2*runif(n); # regressor 
# # 	(i) Homoscedastic erros 
# 	eps <- rnorm(n);	
# #   Heteroscedastic errors 
# 	eps <- rnorm(n)*sigma(x1);  # (ii)
	eps <- (rexp(n)-1)*sigma(x1);		# (iii)
#  	eps <- rexp(n)*sigma(x1);		# (iv)
	Y <- beta1*x1 + eps;		
	
	fit.lm <- lm(Y~x1); 
	
# 	Standard 	confidence intervals 
	temp <- confint(fit.lm)
	CI.stand[i,] <- temp[1,];
	CIs.stand[i,] <- temp[2,];
	
# 	White's estimator of standard deviation 
	se.white <- sqrt(diag(sandwich(fit.lm)));
	
	CI.white[i,] <- fit.lm$coefficients[1] + c(-1,1)*kvant*se.white[1]; 
	CIs.white[i,] <- fit.lm$coefficients[2] + c(-1,1)*kvant*se.white[2]; 
	
# 	White's HC3 estimator of standard deviation 
	se.white.HC3 <- sqrt(diag(vcovHC(fit.lm, type = "HC3")));
	
	CI.HC3[i,] <- fit.lm$coefficients[1] + c(-1,1)*kvant*se.white.HC3[1]; 
	CIs.HC3[i,] <- fit.lm$coefficients[2] + c(-1,1)*kvant*se.white.HC3[2]; 
}




# # CI Standard
# print("Standard confidence intervals:")
b0.stand <- assess.CI(CI.stand, theta=0); 
b1.stand <- assess.CI(CIs.stand, theta=1); 

# # CI White 
# print("White confidence intervals:")
b0.white <- assess.CI(CI.white, theta=0); 
b1.white <- assess.CI(CIs.white, theta=1); 

# # CI White HC3 
# print("White HC3 confidence intervals:")
b0.white.HC3 <- assess.CI(CI.HC3, theta=0); 
b1.white.HC3 <- assess.CI(CIs.HC3, theta=1); 


# # Interecept 
print("CI for intercept:")
round(rbind(b0.stand, b0.white, b0.white.HC3),3);

# # Slope 
print("CI for slope:")
round(rbind(b1.stand, b1.white, b1.white.HC3),3);


# Task: Compare the coverages and lengths of the confidence intervals 
# for the following conditional distribution of errors (given X=x):
# (i) N(0,1)
# (ii) exp(x) * N(0,1)
# (iii) exp(x) * (Exp(1)-1)
# (iv) exp(x) * Exp(1)
# Try to explain the results.  

# # Results for (iii) exp(x) * (Exp(1)-1)
# Number of MC samples 10000

# Results for intercept estimation 
# 
#                      n = 100                   |         n = 500
#              coverage mean(length) sd(length)  | coverage mean(length) sd(length)
# b0.stand        0.991     2.866     0.611      |    0.991    1.27      0.126
# b0.white        0.938     1.960     0.491      |    0.947    0.88      0.110
# b0.white.HC3    0.942     2.018     0.509      |    0.948    0.89      0.111

# Results for slope estimation 
# 
# rbind(b1.stand, b1.white, b1.white.HC3)        |          n = 500
#              coverage mean(length) sd(length)  |  coverage mean(length) sd(length)
# b1.stand        0.896     2.483     0.491      |     0.911     1.105     0.103
# b1.white        0.914     2.817     0.766      |     0.947     1.289     0.175
# b1.white.HC3    0.920     2.904     0.795      |     0.948     1.297     0.176