Prinstats homework 1

Here is the link to get the code from github https://github.com/chamyields/mastats/tree/main/pinstats_2020

Questions

Consider the paper of David Spiegelhalter (available on Ufora). In this paper the author examines the age-specific risk of dying (with or from) Covid-19 in the UK. Table 2 in the paper, gives for several age categories the numbers of reported Covid-19 deaths, numbers of non-covid related deaths, population sizes, among other numbers.

1. Question

Assign to each age interval a numerical value given by the midpoint (e.g. for 0-4, the midpoint is (0+4)/2 =2) and let \(p_i\) denote the corresponding population size in the UK. Assume that the number of non-covid deaths \((Y^N )\) within an age class with midpoint \(a_i\) \((i = 1,\ldots , 11)\) can be described by a Poisson distribution with mean given by

\[ \text{E}(Y^N|\text{age}=a_i) = p_i\exp(\beta_0^N+\beta_1^Na_i) \]

(the superscript refers to non-covid). Also for the covid-related deaths \((Y^C)\) we assume a Poisson distribution with

\[ \text{E}(Y^C|\text{age}=a_i) = p_i\exp(\beta_0^C+\gamma\beta_1^Na_i) \]

How do you interpret the \(\gamma\) parameter?

Solution

Besides the first and the last age intervals, the difference between the mid point from one age group to the other is 10 years.

• Lets take a unit increase in age as 10 years. Then \[ \beta_1^N = \log \text{E}(Y^N|\text{age} = a_i+1) - \log \text{E}(Y^N|\text{age} = a_i) \] and \[ \gamma\beta_1^N = \log \text{E}(Y^C|\text{age} = a_i+1) - \log \text{E}(Y^C|\text{age} = a_i) \]
• \(\beta_1^N\) is the expected change in the number of non-covid-19 deaths per unit increase in age in the logarithmic scale

• \(\gamma \beta_1^N\) is the expected change in the number of covid-19 deaths per unit increase in age in the logarithmic scale

• \(\gamma = \gamma\beta_1^N/\beta_1^N\) is ratio of the expected change in the number of covid-19 deaths on the expected change in the number of non-covid 19 deaths for a unit increase in age in the logarithmic scale
• \(\gamma>0\) implies for a unit increase in age the expected change in the number of covid-19 deaths is greater than that of non-covid 19 deaths.
  
• \(\gamma<0\) implies for a unit increase in age, the expected change in the number of covid-19 deaths is less than that of non-covid 19 deaths.

• \(\gamma = 0\) implies for a unit increase in age, the expected change in the number of covid-19 deaths is equal to that of non-covid 19 deaths. Hence dying from covid 19 compared to dying from a non-covid 19 cause is independent of age.

2. Question

Given the models from the previous question, and assuming that all num- bers of deaths (both covid and non-covid) are independently distributed. Construct the log-likelihood function and make a graph of the log-likelihood as a function of the \(\gamma\) parameter. You may use the following estimates for the other parameters:

\[ \beta_0^N = −12.06 ~~ \beta_0^C = −14.15 ~~ \beta_1^N = 0.09836 \]

Solution

We have been told that \(Y^N\) and \(Y^C\) follow a Poisson distribution, recall that if a random variable \(K\) follows a Poission distribution, it has the follow mass function

\[ P_\lambda(K=k) = f(k|\lambda) = \frac{\lambda^{k}e^{-\lambda}}{k!} \] Let’s say we have observed \(n\) values from \(n\) random variables \(K_1,\ldots,K_n\) that are:

• all identically distributed (i.e. they all follow a Poisson distribution with parameter \(\lambda_i\), well not exactly identical since their mean \(\lambda\)s are different)
• independent, i.e, \(P(K_1,\ldots,K_n) = \prod_i P(K_i)\)

The likelihood is then defined as a function of \(\mathbf \Lambda = (\lambda_1,\ldots,\lambda_n)\) given the data, i.e., \[ L(\mathbf \Lambda|K_1=k_1,\ldots,K_n=k_n) = \prod_i P_{\lambda_i}(K_i=k_i) = \prod_if(k_i|\lambda_i) = \prod_i\frac{\lambda_i^{k_i}e^{-\lambda_i}}{k_i!} \] The log likelihood which is simpler transformation to work with is given by \[ l(\Lambda|K_1=k_1,\ldots,K_n=k_n) = \log \prod_i P_{\lambda_i}(K_i=k_i) = \sum_i\log\frac{\lambda_i^{k_i}e^{-\lambda_i}}{k_i!} \]

\[ l(\Lambda|K_1=k_1,\ldots,K_n=k_n) = \sum_i \big(k_i \log\lambda_i - \lambda_i - \log k_i!\big) \]

\[ l(\Lambda|K_1=k_1,\ldots,K_n=k_n) = \sum_i k_i \log\lambda_i - \sum_i \lambda_i - \sum_i \log k_i! \]

For the problem at hand we have:

• \(\lambda_i^N = \text{E}(Y^N|\text{age}=a_i)\) for non-covid 19 deads
• \(\lambda_i^C (\gamma) = \text{E}(Y^C|\text{age}=a_i)\) for covid 19 deads
• our parameter of interest is \(\gamma\).
• samples size \(n_N=11\) for the covid-19 deads
• samples size \(n_C=11\) for the non-covid-19 deads

\[ l(\gamma|Y^C_1=Y^C_1,\ldots,Y^C_n=Y^C_n,Y^N_1=y^N_1,\ldots,Y^N_n=y^N_n) = \log \Pi_if(Y^C|\lambda_i^C(\gamma))f(y_i^N|\lambda_i^N) \] \[ l(\gamma|Y^C_1=Y^C_1,\ldots,Y^C_n=Y^C_n,Y^N_1=y^N_1,\ldots,Y^N_n=y^N_n) = \sum_{i=1}^{11}\log f(Y^C|\lambda_i^C(\gamma))+\sum_{i=1}^{11}\log f(y_i^N|\lambda_i^N) \] \[ l(\gamma|Y^C_1=Y^C_1,\ldots,Y^C_n=Y^C_n,Y^N_1=y^N_1,\ldots,Y^N_n=y^N_n) = \sum_{i=1}^{11}\big(Y^C\log \lambda_i^C(\gamma) + y_i^N\log \lambda_i^N\big) -\sum_{i=1}^{11}\big(\lambda_i^C(\gamma)+\lambda_i^N\big)-\sum_{i=1}^{11}\log (Y^Cy_i^N) \] We can simply write our log likelihood function as

\[ l(\gamma) = \sum_{i=1}^{11}\big(Y^C\log \lambda_i^C(\gamma) + y_i^N\log \lambda_i^N\big) -\sum_{i=1}^{11}\big(\lambda_i^C(\gamma)+\lambda_i^N\big)-\sum_{i=1}^{11}\log (Y^CY^N) \] So we have our log likelihood function lets code it in R.

library(dplyr,
        verbose=FALSE)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(tidyverse,
        verbose = FALSE)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✔ ggplot2 3.2.1     ✔ purrr   0.3.3
✔ tibble  2.1.3     ✔ stringr 1.4.0
✔ tidyr   1.0.2     ✔ forcats 0.5.0
✔ readr   1.3.1     
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

age_mid_point = function(x,y){
  (x+y)/2
}

beta_0_C_hat = -14.15
beta_1_N_hat = 0.09836
beta_0_N_hat = -12.06

lb = c(0,5,15,25,35,45,55,65,75,85)
ub = c(4,14,24,34,44,54,64,74,84,89)

age = c(age_mid_point(lb,ub),90)

Y_C = c(3,3,33,128,369,1283,3476,7319,16043,10160,10790)

population = c(3535430,7159102,6988755,
               7998302,7460856,8142528,
               7019590,5906928,3476922,
               918437,528959)

Y_N = c(848,138,470,1141,2577,6577,13565,27184,47729,30703,30703)

data = data.frame(age,Y_N,Y_C,population)


beta_0_C = beta_0_C_hat
beta_0_N = beta_0_N_hat
beta_1_N = beta_1_N_hat


log_l = function(gamma_hat){

  lambda_N = data$population*exp(beta_0_N+beta_1_N*data$age)
  lambda_C = data$population*exp(beta_0_C+gamma_hat*beta_1_N*data$age)
  
  sum_1 = data$Y_C*log(lambda_C)+data$Y_N*log(lambda_N)
  sum_2 = lambda_C+lambda_N
  sum_3 = log(data$Y_C*data$Y_N)
  
  sum_all= sum(sum_1-sum_2-sum_3)
  return(-sum_all)
}
data

log_ls = c()
gammas = seq(1,1.2,length.out = 1000)
for(gamma in gammas){
  log_ls  = c(log_l(gamma),log_ls)
}

log_ls = log_ls/max(abs(log_ls))

plot(x=gammas,y=-log_ls,type="l")
abline(v=gammas[-log_ls==max(-log_ls)])

paste0('optimal gamma: ',round(gammas[-log_ls == max(-log_ls)],3))

[1] "optimal gamma: 1.072"

3.Question

Use the optim (or similar) function in R for finding the maxi- mum likelihood estimate of \(\gamma\).

Solution

Lets look for \(\gamma\) with optim

# initial value of gamma
par = 0

opt = optim(par=par,
      fn= log_l)

one-dimensional optimization by Nelder-Mead is unreliable:
use "Brent" or optimize() directly

opt

$par
[1] 1.127734

$value
[1] -1888745

$counts
function gradient 
      32       NA 

$convergence
[1] 0

$message
NULL

Can you compare the results of the manual computation with that of optim. They seem to be different? Please check if there are any mistakes in my computations

plot(x=gammas,y=-log_ls,type="l")
abline(v=gammas[-log_ls==max(-log_ls)],col='blue')

abline(v=opt$par,col='red')

• Note that our log_l function returns the negative log likelihood instead because optim is a minimisation function.

• I checked the paper for the maximum likelihood estimator of \(\gamma\), they also have \((1.128)\), the value from the optim function. So emm…. let’s check the optimise function ;)

optimize(log_l,c(-2,2))

$minimum
[1] 1.127681

$objective
[1] -1888745

4. Question

with the information in Table 2 and with the models and their parameter estimates of the previous questions, estimate the following two probabilities: \[ \text{P} (\text{Age} = 19.5 | \text{ died of covid}) \text{ and } \text{P} (\text{Age} = 79.5 | \text{ died of covid}) \]

Solution

Bayes theorem: for two random variables \(X\) and \(Y\), \[ P(X|Y) = \frac{P(Y|X)P(X)}{P(Y)} \]

Hence

\[ \text{P} (\text{Age} = 19.5 | \text{ died of covid}) = \frac{\text{P} ( \text{ died of covid} | \text{Age} = 19.5 ) \text{P} (\text{Age} = 19.5)}{\text{P}( \text{ died of covid})} \]

• further, use these values and plug into the formula above.

\[ \text{P} ( \text{ died of covid} | \text{Age} = 19.5 ) = \frac{\text{E}(Y^C|\text{Age} = 19.5)}{\text{population size with mid point at age 19.5}} \] \[ \text{P} (\text{Age} = 19.5) = \frac{\text{Population size with mid point at 19.5}}{\text{Total population size across all age groups}} \]

\[ \text{P} (\text{ died of covid}) = \frac{\text{Total number of people who died from covid across all age groups}}{\text{Total population size across all age groups}} \]

Fill the values and compute this probability, do the same for \(\text{P} (\text{Age} = 79.5 | \text{ died of covid})\)

Question 5.

For this question you may assume that the numerical values for (the estimates) \(\beta_0^C\), \(\beta^N_0\) and \(\beta_1^N\) that were given in question 2, are the true population parameter values.Set up a simulation study in which you simulate data (as in Table 2) for the covid and non-covid related deaths, and calculate the maximum likelihood estimate of \(\gamma\) for each simulated dataset.

• Plot a histogram of the estimates of \(\gamma\). How do you interpret this histogram?
• Assess whether the estimates agree with a normal distribution.

Solution

n=11
gamma_hat = 1.2

lambda_N = data$population*exp(beta_0_N+beta_1_N*data$age)
lambda_C = data$population*exp(beta_0_C+gamma_hat*beta_1_N*data$age)

max_lik_est = c()
num_simulations = 1000

for(k in 1:num_simulations){
  
  # set.seed(k)
  
  Y_C = rpois(n,lambda_C)
  Y_N = rpois(n,lambda_N)
  data = data.frame(age,Y_N,Y_C,population)
  
  par = 1
  
  opt = optimize(log_l,c(-1.5,1.5));
  
  max_lik_est[k] = opt$minimum
}
m<-mean(max_lik_est)
std<-sqrt(var(max_lik_est))

hist(max_lik_est, 
     col = 'skyblue3',
     breaks = 20,
     prob=TRUE)
curve(dnorm(x, mean=m, sd=std), 
      col="red", 
      lwd=3, 
      add=TRUE, 
      yaxt="n")

print('')

Assessing normality

qqnorm(max_lik_est, pch = 1, frame = FALSE)
qqline(max_lik_est, col = "steelblue", lwd = 2)

Ok it seems like there are two sets of maximum likelihood estimates:One that occurs at the upper bound of the interval of possible \(\gamma\) values and the others. The upper bound of the interval is taken as the maximum likelihood estimate only when it can not be computed. We will take the upper bound values out of the distribution and reassess.

### taking out max_lik that are greater 1.4
sub_max_lik = max_lik_est[max_lik_est<1.4]

m<-mean(sub_max_lik)
std<-sqrt(var(sub_max_lik))

hist(sub_max_lik, 
     col= 'skyblue3',
     breaks = 20,
     prob=TRUE)

curve(dnorm(x, mean=m, sd=std), 
      col="red", lwd=3, add=TRUE, yaxt="n")
abline(v=mean(sub_max_lik),col="pink",lwd=3)
abline(v=gamma_hat,col="green",lwd=3)

qqnorm(sub_max_lik, pch = 1, frame = FALSE)
qqline(sub_max_lik, col = "steelblue", lwd = 2)

Now the distribution of mles (maximum likelihood estimates) seem to be normal.

Question 6

For this question you may again assume that the numerical values for \(\beta_0^C\), \(\beta^N_0\) and \(\beta_1^N\) that were given in question 2, are the true population parameter values. We now consider the log-likelihood function of question 2 as a function of \(\gamma\), say \(l(\gamma)\). In likelihood theory there is a result that says that in large samples, the following approximately holds:

\[ 2(l(\hat \gamma) − l(\gamma)) \stackrel{d}{=}\chi_1^2 \]

Upon using this result, construct a 95% confidence interval for the parameter \(\gamma\) and report this interval for the data from Table 2 (i.e. for the fitted model from question 3). How do you interpret this result and what is your final conclusion of your data analysis?

Solution

Note that \(2(l(\hat \gamma) − l(\gamma)) \stackrel{d}{=}\chi_1^2\) gives rise to the well known likelihood-ratio (LR) test.

In the LR test of the null hypothesis \[ H_0:\gamma = \gamma_0 \] versus the two-sided alternative \[ H_1:\gamma \neq \gamma_0 \] we would reject \(H_0\) at the \(\alpha\)-level if the LR statistic exceeds the \(100(1-\alpha/2)\)th or is below \(\alpha/2\) percentile of the \(\chi_1^2\) distribution. These values for \(\alpha=0.05\)-level are:

alpha = 0.05
lb = alpha/2
ub = 1-(alpha/2)
# quantiles
qchisq(c(lb,ub),1)

For our problem at with \(\hat \gamma = 1.28\) and let \(\gamma_0 = 0\)

# Note that the log_l function returns the negative log_l and the test statistics (LR)
# here is 
-2*(log_l(1.28)-log_l(0))

This value is is above the upper bound and hence \(\hat \gamma = 1.28\) is significantly non-zero.

We can construct the \(\alpha\) confidence interval as \[ \chi_{1,(\alpha/2)}^2\leq2(l(\hat \gamma) − l(\gamma)) \leq \chi_{1,(1-\alpha/2)}^2 \]

or \[ \frac{\chi_{1,(\alpha/2)}^2}{2}\leq (l(\hat \gamma) − l(\gamma)) \leq \frac{\chi_{1,(1-\alpha/2)}^2}{2} \]

which

qchisq(c(lb,ub),1)/2

Intepret this interval please!