###
###  MATHEMATICAL STATISTICS 3
###  Sandwich estimator of variance 
###

rm(list=ls());

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

load("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"))

points(DATA$concentration, fit$fitted.values, col="blue", pch=8)

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 
round((1/n)*solve(Gammahat),6); 

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

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

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



                                         
# # # A small simulation experiment 
nopak <- 1000;
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 <- 1000;  # 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);
# 	Y <- rbinom(n, size=200, prob=exp(X%*%beta)/200);
# 	Y <- rgeom(n, 1/(1+exp(X%*%beta)));
	
	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)