gamma_fits.Rmd

---
title: "Fitting the Gamma model to Blood pressure data"
author: "David Steinsaltz"
header-includes: \usepackage{commath}
output:
  pdf_document: default
  html_document: default
---

## Distribution of BP measurements
```{r setup, include=FALSE}
require(knitr)
knitr::opts_knit$set(root.dir = '../..')
require(rstan)
require(tidyverse)
require(magrittr)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
```

### Extract intervals for the digits
Suppose the fractions of digits 0,2,4,6,8 are $b_0,b_2,b_4,b_6,b_8$.
Letting $B_0=0$ and $B_k=10\sum_{j=0}^{k-1}b_{2j}$ for $k=1,\dots,5$,
we want to choose a positive $a$ and place breaks at $-a+B_k$, so that measurements between $-a+B_k$ and $-a+B_{k+1}$ modulo
10 are assigned the final digit $2k$, for $k=0,\dots,4$.
We choose $a$ to minimise the total distance of the intervals from the rounded value:
$$
  \sum_{k=0}^4 \int_{-a+B_k}^{-a+B_{k+1}} \bigl| x-2k\bigr|\mathrm{d} x=\frac12\sum_{k=0}^4 \left(-a+B_k-2k\right)^2 + \left(-a+B_{k+1}-2k\right)^2,
$$
as long as $2k$ is in the appropriate interval. This is minimized at
$$
  a= \frac{1}{5}\left(B_1+B_2+B_3+B_4 - 15\right)=\sum_{j=0}^3 (8-2j) b_{2j} \, - 3.
$$
## Load data
```{r load_data, include=FALSE}
w<-getwd()
L1 <- paste0('/Users/davidsteinsaltz/Library/CloudStorage/OneDrive-Nexus365/research/demography/demography papers/NHANES BP/Nhanes2021','/Data_cleaned/nhanesA.RData')
S1 <- paste0('/Users/davidsteinsaltz/Library/CloudStorage/OneDrive-Nexus365/research/demography/demography papers/NHANES BP/Nhanes2021','/NHANESDiagnostics/parameters.R')
load(L1)
source(S1)
```
```{r process, include=FALSE}
N=dim(nhanesA)[1]
k=3
BP_type_names <- c('Systolic','Diastolic')
BP_place_names <- c('Home','Clinic')
whichsys=match(c('systolicA','systolicB','systolicC'),names(nhanesA))
whichdias=match(c('diastolicA','diastolicB','diastolicC'),names(nhanesA)) 

whichsyshome=match(c('systolicAhome','systolicBhome','systolicChome'),names(nhanesA))
whichdiashome=match(c('diastolicAhome','diastolicBhome','diastolicChome'),names(nhanesA))
whichBP=c(whichsys,whichdias)
whichBPhome=c(whichsyshome,whichdiashome)
# allsys=c(whichsys,whichsyshome)
# alldias=c(whichdias,whichdiashome)

# Make BP measures into array
sys=data.matrix(nhanesA[,whichsys])
dias=data.matrix(nhanesA[,whichdias])
sysH=data.matrix(nhanesA[,whichsyshome])
diasH=data.matrix(nhanesA[,whichdiashome])
L=length(sys)
gamma_dimnames <- list( c('alpha','theta','beta') , c('Clinic' , 'Home'))
norm_dimnames <- c('m_M','m_Delta', 'sigma2_M', 'sigma2_Delta')
```
```{r imputation, eval=TRUE, tidy=TRUE}
digitbreaks <- function(bp){
  digit_table <- unname(table(unlist(bp)%%10))
  digits=digit_table/sum(digit_table)
  a <- sum(seq(8,0,by = -2)*digits) - 3 # Find the starting point for the first interval that minimises
  list(BP=bp,digit.breaks=c(0,cumsum(digits[-5]))*10-a)
}

# Combine the measures into a list, so they can be processed uniformly
#   First level, Systolic or Diastolic
#   Second level, Home or Clinic
#  Third level, BP and breaks
allBP <- list(Systolic=list(Clinic=digitbreaks(sys),Home=digitbreaks(sysH)),Diastolic=list(Clinic=digitbreaks(dias),Home=digitbreaks(diasH)))

  # Shift by the mean of the cumulative sums; apply transposes, so we transpose back

### Impute fractional parts. 

impute <- function(d, breaks=c(1,3,5,7,9)){
  # Intervals defined relative to the centre
  right_breaks <- breaks-seq(0,8,by=2)
    left_breaks <- c(breaks[5]-10,breaks[1:4]-seq(2,8,by=2))
  if (!all(d%%2==0)){stop('Not all even digits')}
  else{
    if (is.null(dim(d))){ # Not an array
      d2 <- (d%%10)/2+1
      runif(length(d2))*(right_breaks[d2]-left_breaks[d2])+left_breaks[d2]+d
    }
    else{
      apply(d,2,function(dd) impute(dd,breaks))
    }
  }
}

# Input is a list BP=matrix of measures, digit.breaks=break points
imputeBP <- function(bp_with_breaks){
  breaks <- bp_with_breaks$digit.breaks
  d <- bp_with_breaks$BP
  list(BP=impute(d,breaks),digit.breaks=breaks)
}

allBP_imp <- lapply(allBP, function(BPtype) lapply(BPtype,function(BPplace) imputeBP(BPplace)))
```
## Fit the BP distribution parameters
We suppose that each individual has BP measures
$\tilde{y}_{ij}^l$ for $i=1,\dots,n$, $j=1,\dots,k$ (default $k=3$),
and $l=1,2$, which are rounded versions of 
$$
  y_{ij}^l \sim \mathcal{N}\bigl( \mu_i^l,(\tau_i^l)^{-1}\bigr),
$$
where 
\begin{align*}
\mu_i^1&=(M_i+\Delta_i)/2,\\
\mu_i^2&=(M_i-\Delta_i)/2,\\
M_i&\sim \mathcal{N}\bigl(m_M,\sigma^2_M) \text{ and }
  \Delta_i\sim \mathcal{N}\bigl(m_\Delta,\sigma^2_\Delta) \text{ independent,}\\
  \tau_i^l &\sim \mathrm{Gamma}(\alpha^l,\alpha^l/\theta^l ).
\end{align*}
(Note that $\alpha^l$ is the usual shape parameter,
while $\theta^l$ is the expectation.)

We wish to estimate the eight parameters 
$$
(m_M,m_\Delta,\sigma^2_M,\sigma^2_\Delta,\alpha^1,\theta^1,\alpha^2,\theta^2)
$$
We begin by assuming $y_{ij}^l$ observed directly. We estimate
by maximising the partial likelihood on the observations
\begin{align*}
  \bar{y}_{i+}&:= \frac{1}{2k} \sum_{j=1}^k y_{ij}^1 + y_{ij}^2,\\
  \bar{y}_{i-}&:= \frac{1}{2k} \sum_{j=1}^k y_{ij}^1 - y_{ij}^2,\\
  s_i^l&:=  \frac{1}{k-1}\sum_{j=1}^k \Bigl( y_{ij}^l - \frac{1}{k} \sum_{j=1}^k y_{ij}^l \Bigr)^2.
\end{align*}
Note that 
$$
(k-1)s_i^l \tau_i^l =\sum_{j=1}^k \Bigl( z_{ij}^l - \frac{1}{k} \sum_{j=1}^k z_{ij}^l \Bigr)^2.
$$ 
where $z_{ij}^l$ are i.i.d.\ standard normal
is independent of $\tau_i^l$, thus has a chi-squared distribution
with $k-1$ degrees of freedom --- hence $\frac{k-1}{2}\cdot s_i^l\tau_i^l$ is
gamma distributed with parameters $(\frac{k-1}{2},1)$. Since $\frac{\alpha}{\theta}\tau_i^l$ is independent of $s_i^l\tau_i^l$, with $\mathrm{Gamma}(\alpha,1)$ distribution, we see that $\frac{\theta(k-1)}{2\alpha}s_i^l$ is the ratio of two independent gamma random variables, hence has beta-prime distribution with parameters $\left(\frac{k-1}{2}, \alpha \right)$, so log partial likelihood
$$
  \ell_{\operatorname{Beta}}(\alpha,\theta;s^l_\cdot)=n\alpha \log\frac{\alpha}{\theta}+n\log\Gamma\left(\alpha+\frac{k-1}{2}\right)-n\log\Gamma(\alpha)
  + \frac{k-1}{2} \sum_{i=1}^n \log s_i^l -\left(\alpha+\frac{k-1}{2}\right) \sum_{i=1}^n \log \left(s_i^l+\frac\alpha\theta\right).
$$

The partial Fisher Information has entries
\begin{align*}
 -\frac{\partial^2 \ell}{\partial \alpha^2} &=
   n\psi_1\left(\alpha\right) - n\psi_1\left(\alpha+\frac{k-1}{2}\right)
  - \frac{n}{\alpha} +\sum_{i=1}^n \frac{2\theta s_i^l + \alpha-(k-1)/2}{(\theta s_i^l + \alpha)^2}\\
-\frac{\partial^2 \ell}{\partial \theta^2} &=
   -\frac{n \alpha}{\theta^2} +\frac{\alpha}{\theta^2}\left(\alpha+\frac{k-1}{2}\right)\sum_{i=1}^n \frac{2\theta s_i^l + \alpha}{(\theta s_i^l + \alpha)^2}\\
-\frac{\partial^2 \ell}{\partial \theta\partial\alpha} &= \frac{n}{\theta}-
   \frac1\theta \sum_{i=1}^n \frac{\alpha^2+2\alpha\theta s_i^l+\frac{k-1}{2}\theta s_i^l}{(\theta s_i^l + \alpha)^2}.
\end{align*}
where $\psi_1$ is the trigamma function.
```{r log_likelihoods, tidy=TRUE}
  beta_prime_LL=function(alpha,theta,s,k=3){
    k1= (k-1)/2
    n<- length(s)
    alpha * n * log( alpha / theta ) - n * lbeta(k1 , alpha )+ (k1-1)* sum(log(s)) -
      (alpha + k1 )*sum(log(s+alpha/theta))
  }

beta_prime_gradient= function(alpha,theta,s,k=3){
  n<- length(s)
  k1 <- (k-1)/2
  d_a <- -n*log(alpha/theta)  - n*digamma(alpha) + n*digamma(alpha+ k1 ) - (alpha+ k1) * sum(log(s * theta -k1)/(s*theta+alpha))
  d_t <- -n*alpha/theta + (alpha+ k1 ) * alpha/theta * sum( 1 /(theta*s+ alpha))
  c(d_a,d_t)
}

# Fisher Information
beta_prime_FI=function(alpha,theta,s,k=3){
    n<- length(s)
    k1 <- (k-1)/2
    d_aa <- -n/alpha + sum(( 2*theta*s+alpha - k1 )/(theta*s+alpha)^2) -
      n*trigamma(alpha+k1 ) + n * trigamma(alpha)
    d_tt <- -n* alpha /theta^2 + ( alpha + k1 )*alpha/theta^2 * sum((alpha+2*theta*s)/
                                           (theta*s+alpha)^2) 
    d_at <- sum((s^2*theta - s*k1) / (alpha+ theta*s)^2)
    matrix(c(d_aa,d_at,d_at,d_tt), 2, 2)
}

alphatheta <- function(BP){
  k <- dim(BP)[2]
  s <- (k-1)/2*apply(BP,1,var)
  LL <- function(gamma_params){
    -beta_prime_LL(gamma_params[1],gamma_params[2],s,k)
  }
  LL_grad <- function(gamma_params){
    -beta_prime_gradient(gamma_params[1],gamma_params[2],s,k)
  }
  # Using optim because constrOptim doesn't work
  #fit <- suppressWarnings( constrOptim(c(1,1),f = LL, grad = LL_grad, ui = diag(1,nrow = 2,ncol = 2), ci = c(0,0)) ) # Constraint matrix identity, so constrained >0
  fit <- suppressWarnings( optim(par = c(1,1),fn = LL, gr = LL_grad ) )
  fisher_info <- beta_prime_FI(fit$par[1],fit$par[2],s,k)
  list( alpha = fit$par[1], theta = fit$par[2], variance = solve(fisher_info), LogLikelihood = -fit$value)
}
```
Let $(\hat\alpha^l,\hat\beta^l)$ be the maximum partial likelihood estimators. Conditioned on $(\tau_i^l)$ we have
\begin{align*}
  \bar{y}_{i+}&\sim \mathcal{N}\left(m_M, \sigma^2_M + \frac{1}{4k}\left( \frac{1}{\tau_i^1}+\frac{1}{\tau_i^2}\right)\right),\\
  \bar{y}_{i-}&\sim \mathcal{N}\left(m_\Delta,\sigma^2_\Delta + \frac{1}{4k}\left( \frac{1}{\tau_i^1}+\frac{1}{\tau_i^2}\right)\right).
\end{align*}
We would then have MLEs
\begin{align*}
  \hat{m}_M&= \frac{1}{n} \sum_{i=1}^n \bar{y}_{i+},\\
  \hat{m}_\Delta&= \frac{1}{n} \sum_{i=1}^n \bar{y}_{i-},
\end{align*}
which are approximately normally distributed, with means $m_M$ and $m_\Delta$ respectively, and conditional on $\tau_i^l$ standard errors
$$
  \frac{\sigma_M^2}{n}+\frac{1}{4kn^2} \sum_{i=1}^n (\tau_i^1)^{-1} + (\tau_i^2)^{-1} \quad \text{ and } \quad
  \frac{\sigma_\Delta^2}{n}+\frac{1}{4kn^2} \sum_{i=1}^n (\tau_i^1)^{-1} + (\tau_i^2)^{-1},
$$
which we may approximate --- with error on the order of $n^{-3/2}$ --- replaceing the mean of $(\tau_i^l)^{-1}$ by its expected value $\beta^l/(\alpha^l-1)$ to obtain
\begin{align*}
  \mathrm{Var}(\hat{m}_M) &\approx \frac{\sigma_M^2}{n}+\frac{1}{4kn}\left( \frac{\beta^1}{\alpha^1-1}+ \frac{\beta^2}{\alpha^2-1}\right) \\
  \mathrm{Var}(\hat{m}_\Delta) &\approx \frac{\sigma_\Delta^2}{n}+\frac{1}{4kn}\left( \frac{\beta^1}{\alpha^1-1}+ \frac{\beta^2}{\alpha^2-1}\right) 
\end{align*}
Finally, conditioned on the $\tau_i^l$ we have that the random variables $\bar{y}_{i+}$ are normal with variance
$$
  \sigma_M^2+\frac{1}{4k}\left((\tau_i^1)^{-1} + (\tau_i^1)^{-1} \right),
$$
so the unconditional variance is the expected value, or
$$
  \sigma_M^2+\frac{1}{4k}\left(\frac{\beta^1}{\alpha^1-1}+ \frac{\beta^2}{\alpha^2-1} \right).
$$
This yields the estimators
\begin{align*}
  \hat\sigma_M^2 &=\frac{1}{n-1}\sum_{i=1}^n\left(\bar{y}_{i+}-n^{-1}\sum_{i=1}^n y_{i+}\right)^2 - \frac{1}{4k}\left(\frac{\hat\beta^1}{\hat\alpha^1-1}+ \frac{\hat\beta^2}{\hat\alpha^2-1} \right),\\
  \hat\sigma_\Delta^2 &=\frac{1}{n-1}\sum_{i=1}^n\left(\bar{y}_{i-}-n^{-1}\sum_{i=1}^n y_{i-}\right)^2 - \frac{1}{4k}\left(\frac{\hat\beta^1}{\hat\alpha^1-1}+ \frac{\hat\beta^2}{\hat\alpha^2-1} \right).
\end{align*}
Using the delta method, and the fact that  we see that the variance of $\hat\beta/(\hat\alpha-1)$ is approximately
$$
  \frac{\sigma_\beta^2}{(\hat\alpha-1)^2} + \frac{\hat\beta^2\sigma_\alpha^2}{(\hat\alpha-1)^4},
$$
where $\sigma_\alpha$ and $\sigma_\beta$ are the standard errors for $\hat\alpha$ and $\hat\beta$ respectively, so the standard errors for $\hat\sigma_M^2$ and $\hat\sigma_\Delta^2$ are approximately
\begin{align*}
  \operatorname{SE}\left(\hat\sigma_M^2\right)&\approx \frac{1}{2k}\Bigl(\frac{8k^2\hat\sigma_M^2}{n} + \frac{\sigma_\beta^2}{(\hat\alpha^1-1)^2} + \frac{(\hat\beta^1)^2\sigma_\alpha^2}{(\hat\alpha^1-1)^4} + \frac{\sigma_\beta^2}{(\hat\alpha^2-1)^2} + \frac{(\hat\beta^2)^2\sigma_\alpha^2}{(\hat\alpha^2-1)^4} \Bigr)^{1/2},\\
  \operatorname{SE}\left(\hat\sigma_\Delta^2\right)&\approx \frac{1}{2k}\Bigl(\frac{8k^2\hat\sigma_\Delta^2}{n} + \frac{\sigma_\beta^2}{(\hat\alpha^1-1)^2} + \frac{(\hat\beta^1)^2\sigma_\alpha^2}{(\hat\alpha^1-1)^4} + \frac{\sigma_\beta^2}{(\hat\alpha^2-1)^2} + \frac{(\hat\beta^2)^2\sigma_\alpha^2}{(\hat\alpha^2-1)^4} \Bigr)^{1/2}
\end{align*}
```{r parameter_fitting, tidy=TRUE, cache=TRUE,echo=FALSE, eval=TRUE}
# Input a pair of BP lists (clinical and home), and output 8 parameters, with variances for each
all_BP_parameters <- function(measures){
  stopifnot(setequal(names(measures) , c('Clinic','Home')),all.equal(dim(measures$Clinic),dim(measures$Home)))
  k <- dim(measures$Clinic$BP)[2]
  n <- dim(measures$Clinic$BP)[1]
  s <- lapply(measures,function(m) (k-1)/2*apply(m$BP,1,var))
  ybar <- vapply(measures,function(m) apply(m$BP,1,mean),FUN.VALUE = rep(0,n))
  ybar_plus <- (ybar[,'Home']+ybar[,'Clinic'])/2
  ybar_minus <- (ybar[,'Home']-ybar[,'Clinic'])/2
  m_M <- mean(ybar_plus)
  m_Delta <- mean(ybar_minus)
  gamma_fits <- lapply(measures,function(m) alphatheta(m$BP))
  variance_correction <- sum(sapply(gamma_fits,function(GF) GF$alpha/GF$theta/(GF$alpha - 1)))/4/k
  vv_correction <- sum(sapply(gamma_fits,function(GF) diag(GF$variance)*c(GF$alpha^2/GF$theta^2/(GF$alpha-1)^4,1/(GF$alpha-1)^2))) # correction to the variance of the variance
  sigma2_M <- var(ybar_plus) - variance_correction
  sigma2_Delta <- var(ybar_minus) - variance_correction
  list(Gamma=gamma_fits,Normal=list(M=list(m=m_M,m.std.error=sd(ybar_plus)/sqrt(n),sigma2 = sigma2_M, sigma2.std.error = sqrt(8*k*k*sigma2_M/n+vv_correction)),Delta=list(m=m_Delta,m.std.error=sd(ybar_minus)/sqrt(n),sigma2 = sigma2_Delta, sigma2.std.error = sqrt(8*k*k*sigma2_Delta/n+vv_correction))))
}

# Apply separately to Systolic and Diastolic
BP_parameters <- lapply(allBP_imp,function(bp) all_BP_parameters(bp))
save(BP_parameters, file = 'BP_parameters.RData')

## Define functions to convert parameters between the format where 
##   Sys and Dias are the top level, and where normal and gamma are separate
##   convert_param1 takes list (systolic, diastolic), outputs list (normal,gamma)
##   convert_param2 does the opposite
convert_param1 <- function(all_params){
  norm_params <- lapply(all_params, function(params) {
    with(params$Normal, c(m_M = M$m , m_Delta = Delta$m, sigma2_M = M$sigma2, sigma2_Delta = Delta$sigma2))
  })
  gamma_params <- lapply(all_params, function(params) {
    GP <- array( 0, dim = c(3,2) , dimnames = gamma_dimnames)
  GP[ 'alpha', ] <- c(params$Gamma$Clinic$alpha, params$Gamma$Home$alpha)
  GP[ 'theta', ] <- c(params$Gamma$Clinic$theta, params$Gamma$Home$theta)
  GP[ 'beta', ] <- GP[ 'alpha', ] / GP[ 'theta', ]
  GP
  })
  list(Normal = norm_params, Gamma = gamma_params )
}
```
## Test the effect of imputation

```{r test_imputation, eval=FALSE}
number_imp <- 100 # This is the number of different imputations we'll simulate
alphas <- rep(0,number_imp)
thetas <- rep(0,number_imp)
d <- allBP$Systolic$Clinic
for (i in seq_len(number_imp)){
  d_imp <- imputeBP(d)$BP
  ab <- alphatheta(d_imp)
  alphas[i] <- ab$alpha
  thetas[i] <- ab$theta
}
```
## Test whether parameters are being fit correctly
Simulate bootstrap data sets. Find the average parameter estimate, and compare to the "true" parameters.
Also compare the average estimated SE to the "true" SE (which is the SD of the estimates).
```{r test_parameters, eval=TRUE, echo=FALSE, cache = TRUE, results= 'asis'}
k <- 3
num_reps <- 100

BP_simulation <- function(norm_params,gamma_params,num_indiv, num_sim=1,k=3){
  # norm_params is a vector with named entries 'm_M','m_Delta','sigma_M','sigma_Delta',
  # gamma_params is a matrix with rows 'alpha' and 'beta', columns 'Home" and 'Clinic"
  # n is number of simulations
   stopifnot(all(is.element(c('alpha','beta'),dimnames(gamma_params)[[1]])),
              all(is.element(c('Clinic','Home'),dimnames(gamma_params)[[2]])),
              setequal(names(norm_params),norm_dimnames)
    )
  lapply(seq_len(num_sim),function(i){
    M <- rnorm(num_indiv,mean = norm_params['m_M'],sd = sqrt(norm_params['sigma2_M']))
    Delta <- rnorm(num_indiv,mean = norm_params['m_Delta'], sd = sqrt(norm_params['sigma2_Delta']))
    muH <- (M+Delta)
    muC <- (M-Delta)
    tauH <- rgamma(num_indiv,shape = gamma_params['alpha','Home'], rate = gamma_params['beta','Home'])
    tauC <- rgamma(num_indiv,shape = gamma_params['alpha','Clinic'], rate = gamma_params['beta','Clinic'])
    all_individual_params=cbind(mu.Home = muH, mu.Clinic = muC, sigma.Home = tauH^(-.5), sigma.Clinic = tauC^(-.5))
    list(Home = list(BP = t(apply(all_individual_params, 1, function(pars) rnorm(k,pars['mu.Home'],pars['sigma.Home'])))),
             Clinic = list(BP= t(apply(all_individual_params, 1, function(pars) rnorm(k,pars['mu.Clinic'],pars['sigma.Clinic'])))))
  })
}

# Simulation that inputs two sets of parameters, and outputs pairs (systolic, diastolic)
BP_simulation2 <- function(norm_params,gamma_params, num_indiv, num_sim=1,k=3){
  stopifnot(setequal(names(norm_params) , BP_type_names), 
            setequal(names(gamma_params) , BP_type_names))
  bp <- lapply(BP_type_names, function(BPtype) BP_simulation(norm_params[[BPtype]],gamma_params[[BPtype]],num_indiv, num_sim,k))
  names(bp) <- BP_type_names
  lapply(seq_len(num_sim),function(i) lapply(BP_type_names, function(BP_type) bp[[BP_type]][[i]]) %>% 
           setNames(BP_type_names))
}

round_bp <- function(bp,rounding_boundaries=c(-1, 1,3,5,7), real_digits=c(0,2,4,6,8)){
  stopifnot(length(rounding_boundaries)==length(real_digits))
  if (is.array(bp)){apply(bp,2,round_bp)}
  else{
    bp_shifted <- bp - rounding_boundaries[1] # Shift the lowest boundary up to 0
    bp_int <- bp_shifted %% 10
    bp10 <- bp_shifted - bp_int
    sapply(seq_along(bp),function(i) {
      real_digits[sum(bp_int[i]>= (rounding_boundaries-rounding_boundaries[1]) ) ]} ) +bp10
  }
}

# First simulate BPs, then round, according to a different set of boundaries for each place
BP_simulation_rounded <- function(norm_params,gamma_params,num_indiv, num_sim=1,
                                  k=3,rounding_boundaries=list(Home=c(-1, 1,3,5,7), Clinic= c(-1,1, 3, 5,7)),
                                  real_digits=c(0,2,4,6,8)){
  bp <- BP_simulation(norm_params,gamma_params,num_indiv, num_sim,k)
  lapply(bp,function(simulated_dataset) lapply(BP_place_names, function(BPplace) list(BP=round_bp(simulated_dataset[[BPplace]]$BP,rounding_boundaries[[BPplace]], real_digits),digit.breaks=rounding_boundaries[[BPplace]])) %>%
           setNames(BP_place_names)
         )
}

BP_simulation2_rounded <- function(norm_params,gamma_params, num_indiv,
                                   num_sim=1,k=3, rounding_boundaries=list(Home=c(-1, 1,3,5,7), Clinic= c(-1,1, 3, 5,7)),
                                  real_digits=c(0,2,4,6,8)){
  stopifnot(setequal(names(norm_params) , BP_type_names), 
            setequal(names(gamma_params) , BP_type_names))
  bp <- lapply(BP_type_names, function(BPtype) BP_simulation_rounded(norm_params[[BPtype]],gamma_params[[BPtype]],num_indiv, num_sim,k,rounding_boundaries, real_digits)) %>%
    setNames(BP_type_names)
  lapply(seq_len(num_sim),function(i) lapply(BP_type_names, function(BPtype) bp[[BPtype]][[i]])) %>% setNames(BP_type_names)
}

num_reps <- 100
params <- convert_param1(BP_parameters)
#b= BP_simulation2(norm_params , gamma_params2,num_indiv = N,num_sim = num_reps)


for(which_bp in BP_type_names){
  estimates <- array(0,dim = c(num_reps,8),dimnames = list( NULL , c(norm_dimnames,'alpha_C','theta_C','alpha_H','theta_H')))
  std_errors <- array(0,dim = c(num_reps,8),dimnames = list( NULL , c(norm_dimnames,'alpha_C','theta_C','alpha_H','theta_H')))

  bp <- BP_simulation(params$Normal[[which_bp]],params$Gamma[[which_bp]], num_indiv=N, num_reps)
simulation_params <- lapply(bp, function(bp_sim) all_BP_parameters(bp_sim)) # simulated parameter estimates

for(i in seq_len(num_reps)){
  spe <- simulation_params[[i]]
  estimates[i,'m_M'] <- spe$Normal$M$m
  estimates[i,'m_Delta'] <- spe$Normal$Delta$m
  estimates[i,'sigma2_M'] <- spe$Normal$M$sigma2
  estimates[i,'sigma2_Delta'] <- spe$Normal$Delta$sigma2
  std_errors[i,'m_M'] <- spe$Normal$M$m.std.error
  std_errors[i,'m_Delta'] <- spe$Normal$Delta$m.std.error
  std_errors[i,'sigma2_M'] <- spe$Normal$M$sigma2.std.error
  std_errors[i,'sigma2_Delta'] <- spe$Normal$Delta$sigma2.std.error
  estimates[i, 'alpha_C'] <- spe$Gamma$Clinic$alpha
  estimates[i, 'theta_C'] <- spe$Gamma$Clinic$theta
  estimates[i, 'alpha_H'] <- spe$Gamma$Home$alpha
  estimates[i, 'theta_H'] <- spe$Gamma$Home$theta
  std_errors[i, c('alpha_C','theta_C')] <- sqrt(diag(spe$Gamma$Clinic$variance))
  std_errors[i, c('alpha_H','theta_H')] <- sqrt(diag(spe$Gamma$Home$variance))
}

simulation_mean_estimate <- apply(estimates,2,mean)
simulation_mean_SE <- apply(std_errors,2, mean )
simulation_SE <- apply(estimates,2,sd)
true_params <- simulation_SE  # Just a dummy with the right names
true_params[norm_dimnames] <- params$Normal[[which_bp]][norm_dimnames]
true_params[c('alpha_C','theta_C','alpha_H','theta_H')]<-params$Gamma[[which_bp]][c('alpha','theta'),c('Clinic','Home')]
#Make a matrix for comparing true parameters to average simulation estimates
cat(which_bp , "Parameter estimates",'\n')
print(signif(matrix(c(simulation_mean_estimate,true_params,
                    simulation_mean_estimate/true_params-1), 
                  nrow= 8,ncol = 3,
                  dimnames=list(names(simulation_SE),
                                c("Sim.average",'True','Error'))),3))
cat('\n',which_bp, 'parameter SE','\n')
print(matrix(c(simulation_mean_SE, simulation_SE,
         simulation_mean_SE / simulation_SE -1), nrow= 8,ncol = 3,
           dimnames=list(names(simulation_SE),
                     c("Sim.average",'True SE', 'Error'))))
}
```


## Multiple imputation for the real data
```{r multiple_impute, cache=TRUE, tidy= TRUE}
number_impute <- 10
all_variables_est <- c('m_M','m_Delta','sigma2_M','sigma2_Delta','alpha_C','theta_C','alpha_H','theta_H')
all_variables_SE <- c( setNames(vapply(all_variables_est,function(N) paste0(N,'.SE'),'x'),NULL))
all_variables_cov <- c('Covar_C','Covar_H')
all_variables <- c(all_variables_est, all_variables_SE, all_variables_cov)
impute_results <- lapply(allBP, function(BPtype) 
  setNames(data.frame( array(0,dim = c(number_impute,length(all_variables))) ), all_variables))

suppressWarnings(remove('Delta','Gamma'))

## Make containers for the results, one each for systolic and diastolic
for(i in seq_len(number_impute)){
    allBP_imp <- lapply(allBP, function(BPtype) lapply(BPtype,function(BPplace) imputeBP(BPplace)))
      # quantifying over allBP gives Systolic vs diastolic
      # Next level gives Home vs Clinic
      # Next level has BPs and digit breaks
    for (BPtype in BP_type_names){
      spe <- all_BP_parameters(allBP_imp[[BPtype]]) # simulated parameter estimates
      attach(spe$Gamma)
      attach(spe$Normal)
      gamma_SE_C <- sqrt(diag(Clinic$variance))
      gamma_SE_H <- sqrt(diag(Home$variance))
    
  impute_results[[BPtype]][i,all_variables_est] <- c(M$m, Delta$m, M$sigma2, Delta$sigma2,
                                  Clinic$alpha,Clinic$theta, Home$alpha,Home$theta)
  impute_results[[BPtype]][i,all_variables_SE] <- c(M$m.std.error,Delta$m.std.error, M$sigma2.std.error, Delta$sigma2.std.error,
                                   gamma_SE_C,gamma_SE_H)
  impute_results[[BPtype]][i,all_variables_cov] <- c( Clinic$variance[1,2], Home$variance[1,2] )
    detach(spe$Gamma)
    detach(spe$Normal)
    }
}

total_M_I_results <- lapply(impute_results, function(result) 
{ 
  M_I_estimates = apply(result[,all_variables_est],2,mean)
                M_I_std_error = apply(result[,all_variables_est], 2, sd)
                M_I_covariance = c('Clinic' = cov(result[,'alpha_C'],result[,'theta_C']), 'Home' = cov(result[,'alpha_H'],result[,'theta_H']))
  c(list(estimates= M_I_estimates, std_errors = c(sqrt(apply(result[,all_variables_SE]^2, 2 , mean ) + M_I_std_error^2), apply(result[,all_variables_cov], 2 , mean )+M_I_covariance)))
})

total_M_I_results <- lapply(total_M_I_results, function(result) c(result, list(Correlation = setNames(c( Clinic = result$std_errors['Covar_C']/prod(result$std_errors[c('alpha_C.SE','theta_C.SE')]),
                                                                                           'Home' = result$std_errors['Covar_H']/prod(result$std_errors[c('alpha_H.SE','theta_H.SE')])),c('Corr_C','Corr_H'))))
)
```
# Now compute the combined variance
For a parameter like $\alpha$ we estimate the variance of $\hat\alpha$ by
\newcommand{\E}{\mathbb{E}}
\renewcommand{\P}{\mathbb{P}}
$$
  \mathrm{Var}(\hat\alpha) = \mathbb{E}\bigl[ \mathrm{Var}\left(\hat\alpha\, |\, I\right)\bigr] + \mathrm{Var}\left(\mathbb{E} \left[ \hat\alpha\, |\, I \right]\right).
$$
Here $I$ represents the randomly imputed fractional part. 
We can estimate the first term by averaging the estimated variance (from Fisher Information) over all random imputations.
We estimate the second term by the variance of the $\alpha$ estimates over imputations. Note that this is not quite right, since what we really
want the variance of is $\alpha_0(I)$ --- effectively, the ``true'' parameter consistent with the imputation. This is a plug-in estimate,
as is the Fisher Information estimate of the variance.

## Computing residuals


We define the deviance for an individual $i$ with observations $(Y_i)$
given the hyperparameters $h=(m_M,m_\Delta,\sigma^2_M,\sigma^2_\Delta,\alpha^H,\theta^H,\alpha^C,\theta^C)$
$$
  D= \sum_{i=1}^n \log \mathbb{P}\left\{ \mathbf{Y}_{i}\,|\, \text{hyperparameters}=h\right\}.
$$
\newcommand{\wtb}{\widetilde\mathbf}
Since the $\mathbf{Y}_i$ are independent conditioned on $h$,
\begin{align*}
D&= \sum_{i=1}^n \log \E_h\left[ \P\left\{ \mathbf{Y}_i \, |\, M_i,\Delta_i,\tau_i^{C},\tau_i^H \right\} \right]\\
    &\approx \sum_{i=1}^n \log \frac1R\sum_{r=1}^R \left[ \P\left\{ \mathbf{Y}_i \, |\, M_{i,r},\Delta_{i,r},\tau_{i,r}^{C},\tau_{i,r}^{H} \right\}\right] \frac{\pi_h(M_{i,r},\Delta_{i,r},\tau_{i,r}^{C},\tau_{i,r}^{H} )}{q(M_{i,r},\Delta_{i,r},\tau_{i,r}^{C},\tau_{i,r}^{H} \, | \, h,\, \mathbf{Y}_i)},
\end{align*}
where $(M_{i,r},\Delta_{i,r},\tau_{i,r}^{C},\tau_{i,r}^{H})$ are independent samples from a distribution $q$ that may depend
on $\mathbf{Y}_i$ and $h$, and $\pi_h$ is the true density of those individual parameters given hyperparameters $h$.

<!-- 
We can try estimating this simply by Monte Carlo sampling of the individual parameters.
-->
<!--
We estimate this by importance sampling on the four individual parameters $(M_i,\Delta_i,\tau_i^H, \tau_i^C)$ from an approximate posterior.
We have, conditioned on the observations of sample variances for clinical and home measures $S_{Ci}^2$ and $S_{Hi}^2$,
\begin{align*}
  \tau_i^H &\sim \mathrm{Gamma}\left( \alpha^H+\frac{k-1}{2},\, \beta^H+\frac{k-1}{2} S_{Hi}^2 \right),\\
  \tau_i^C &\sim \mathrm{Gamma}\left( \alpha^H+\frac{k-1}{2},\, \beta^H+\frac{k-1}{2} S_{Hi}^2 \right).
\end{align*}
Then, recalling the definitions of $y_{i+}$ and $y_{i-}$,
conditioned on $\tau_i^H$ and $\tau_i^C$ we have
\begin{align*}
  M_i &= \mathcal{N} \left( (\tau_i^M)^{-1} \left( \frac{m_M}{\sigma_M^2} +\y_{i+}\cdot \frac{4k \tau_i^C \tau_i^H}{\tau_i^C + \tau_i^H}  \right) \, ,  (\tau_i^M)^{-1} \right),\\
  \Delta_i &= \mathcal{N} \left( (\tau_i^Delta)^{-1} \left( \frac{m_\Delta}{\sigma_\Delta^2} +\y_{i+}\cdot \frac{4k \tau_i^C \tau_i^H}{\tau_i^C + \tau_i^H}  \right) \, ,  (\tau_i^\Delta)^{-1} \right),\\
\end{align*}
where
\begin{align*}
  \tau_i^M & = \frac{1}{\sigma_M^2} + \frac{4k \tau^C_i \tau^H_i}}{\tau^C_i + \tau^H_i},\\
  \tau_i^\Delta & = \frac{1}{\sigma_\Delta^2} + \frac{4k \tau^C_i \tau^H_i}}{\tau^C_i + \tau^H_i},\\
\end{align*}
Because of rounding, the observed $S^2_i$ will be too small. We approximate the true $S^2_i$ by adding $\frac13$, the variance
of a uniform random variable on $[-1,1]$.
-->

## Check empirical variances
We simulate new variances from the inferred model, and compare
it with QQ plots. The first thing we do is to compare the empirical
variances (with imputed fractional parts) to the beta-prime distribution.
```{r Variance test}
# Compare to beta-prime distribution
params <- convert_param1(BP_parameters)
dbetaprime <- function(x, A=1, B= BP_parameters$Systolic$Gamma$Clinic$alpha)
  {
    x^(A-1)*(1+x)^(-A-B)/ beta(A,B)
}

all_s2k <- NULL
all_s2k_sim <- NULL
all_s2k_sim_noimp <- NULL
all_SD <-NULL
all_HC <- NULL

oversample <- 1

bp_sim <- BP_simulation2(params$Normal,params$Gamma,
                      num_indiv= N, num_sim = oversample) # Get oversample sets of simulations; 
    #We're going to take 1/oversample of them, to get a more accurate
#    estimate of the distribution, particularly at the upper end

for(which_bp in BP_type_names){
#    s2k <- apply(bp_sim[[which_bp]])
    for(i in seq_len(oversample)){
      allBP_imp <- lapply(allBP, function(BPtype) lapply(BPtype,function(BPplace) imputeBP(BPplace)))
    for ( where_bp in BP_place_names){
      
      s2k <- apply(allBP_imp[[which_bp]][[where_bp]]$BP,1,var) /
                        params$Gamma[[which_bp]]['beta',where_bp]
      s2k_sim <- apply(bp_sim[[i]][[which_bp]][[where_bp]]$BP,1,var) /
                         params$Gamma[[which_bp]]['beta',where_bp]
      all_s2k %<>% c(sort(s2k))
      all_s2k_sim %<>% c(sort(s2k_sim)[seq(from = ceiling(oversample/2), by = oversample, length.out = N)])
      all_SD %<>% c(rep(which_bp,length(s2k)))
      all_HC %<>% c(rep(where_bp,length(s2k)))
}
    ggplot(tibble(Variance = s2k[s2k<4]), aes(x=Variance)) + 
      geom_histogram(alpha=0.2, position="identity", aes(y = ..density..), binwidth = .04) +
      labs(title =  paste(which_bp,where_bp, ': Histogram of Variances' )) +
      stat_function(fun = function(xx) dbetaprime(xx,A =1, 
                          B = params$Gamma[[which_bp]]['alpha',where_bp]))
  }
}

for(which_bp in BP_type_names){
  for(where_bp in BP_place_names){
    all_s2k[all_SD == which_bp & all_HC == where_bp] <- sort(all_s2k[all_SD == which_bp & all_HC == where_bp])
    all_s2k_sim[all_SD == which_bp & all_HC == where_bp] <- sort(all_s2k_sim[all_SD == which_bp & all_HC == where_bp])
  }
}
# Need to sort the variances within simulation groups.

```
Now we compare the distribution of observed variances to the distribution of simulated variances.
```{R Compare variances}
all_s2k.tbl <- subset(tibble(Real = all_s2k, Simulated = all_s2k_sim, 
                      HC = factor(all_HC),
                      DS = factor(all_SD)), 
                      pmax(Real, Simulated) < 12 )

s2k.plot <- ggplot(all_s2k.tbl, aes(x=Real, y=Simulated, colour=interaction(HC,DS),
  group=interaction(HC, DS))) + 
  geom_point() + geom_line() +  stat_function(fun = identity , color = 'black' ) + labs( colour = 'BP type', title = 'Real vs. Simulated variances with imputed fractional parts') +
  scale_color_brewer(palette="Dark2")

s2k.plot
```
## Estimating the mortality parameters.
We divide the individuals up into three races $r_i\in \{\mathrm{B, W, M}\}$
and two sexes $s_i\in \{ \mathrm{M, F}\}$.
For each of these we model the mortality rate at age $t$ as 
$h_i(t)=B_{r_i s_i} \mathrm{e}^{\theta_{r_i s_i} t}$. So we need to estimate

```{r Deviance, echo=FALSE}
# Function to check format of parameters, blood pressure data
#  Must have right 4 parameters for normal and gamma parameters
#   gamma parameters arranged as table, with clinic and home parameters
#   bp must be a list of Clinic and Home values
#   Each element includes a BP item, with a table of BP measures, of the same dimensions

check_parameters <- function(norm_params,gamma_params,bp){
stopifnot(all(is.element(c('alpha','beta'),dimnames(gamma_params)[[1]])),
              all(is.element(c('Clinic','Home'),dimnames(gamma_params)[[2]])),
              setequal(names(norm_params),c('m_M','m_Delta', 'sigma_M', 'sigma_Delta')),
                       is.list(bp), setequal(names(bp),c('Clinic','Home')),
          all(sapply(bp,function(BPplace) is.element('BP',names(BPplace)))),
          dim(bp$Clinic$BP)==dim(bp$Home$BP)
    )
}

# Function to compare variances between true rounded observations bp
#   and variances of simulations from (norm_params, gamma_params)
test_given_parameters <- function(norm_params, gamma_params, bp, num_samples, real_digits = c(0,2,4,6,8)){
  check_parameters(norm_params,gamma_params,bp) # Stop if parameters or bp don't have right format
  stopifnot(all(sapply(bp,function(BPplace) sapply(BPplace$BP %% 10, function(bb) bb %in% real_digits ))))
  num_indiv <- dim(bp$Clinic$BP)[1]
  k <- dim(bp$Clinic$BP)[2]
  S2k <- lapply(bp,function(BPplace) apply(BPplace$BP,1,var)) %>% setNames(BP_place_names)
  mu <- sapply(bp,function(BPplace) apply(BPplace$BP,1,mean))
  MDelta <- t(apply(mu,1, function(mCH) c(mean(mCH),(mCH['Home']-mCH['Clinic'])/2)))
  colnames(MDelta) <- c('M','Delta')
  digit_breaks <-lapply(bp, function(placeBP) placeBP$digit.breaks  ) %>% setNames(BP_place_names)
  sim_data <- BP_simulation_rounded(norm_params, gamma_params, num_indiv, num_samples, k , digit_breaks)
  S2k_simulated <- lapply(BP_place_names, function(placeBP) sapply(sim_data, function(SimD) apply(SimD[[placeBP]]$BP,1,var))) %>%
    setNames(BP_place_names)
  list(real=S2k, sim= S2k_simulated)
}

# dev_resid <- function(norm_params,gamma_params,bp,num_samples,real_digits=c(0,2,4,6,8)){
#   check_parameters(norm_params,gamma_params,bp) # Stop if parameters or bp don't have right format
#   stopifnot(all(sapply(bp,function(BPplace) sapply(BPplace$BP %% 10, function(bb) bb %in% real_digits ))))
#   num_indiv <- dim(bp$Clinic$BP)[1]
#   k <- dim(bp$Clinic$BP)[2]
#   # Make a list of 'Home' and 'Clinic' entries, each of which has an
#   # upper and lower bound for the observed rounded BP.
#   upper_lower <- lapply(bp,function(BPplace) {
#     rounding_boundaries <- BPplace$digit.breaks
#     lapply(list(lower = c(rounding_boundaries[length(rounding_boundaries)]-10, BPplace$digit.breaks[-length(rounding_boundaries)]), 
#                   upper = rounding_boundaries), function(bounds)
#       10* (BPplace$BP %/% 10) + apply(BPplace$BP,2,function(bb) 
#         bounds[match(bb %% 10, real_digits) ] ) ) 
#     }
#   ) # First create a list with lower interval bounds and upper interval bounds;
#   # Then find the interval for each last digit and do the rounding.
#   S2k <- sapply(bp,function(BPplace) apply(BPplace$BP,1,var) + 1/3)*(k-1)/2
#   mu <- sapply(bp,function(BPplace) apply(BPplace$BP,1,mean))
#   MDelta <- t(apply(mu,1, function(mCH) c(mean(mCH),(mCH['Home']-mCH['Clinic'])/2)))
#   colnames(MDelta) <- c('M','Delta')
#   tau <- array(0,dim = c(num_indiv,num_samples,2),dimnames = list(NULL,NULL,BP_place_names))
#   means <- array(0,dim = c(num_indiv,num_samples,4),dimnames = list(NULL,NULL,c('M','Delta','Home','Clinic')))
#   # Simulate tau values for everyone
#   for(place in BP_place_names){
#     tau[,,place] <- rgamma(num_samples*num_indiv,shape= gamma_params['alpha',place]+(k-1)/2, 
#                            rate = gamma_params['beta',place]+S2k[,place])
#   }
#   # Precision for the posterior mean
#     tau_combined <- 4*k*tau[,,'Clinic']*tau[,,'Home']/(tau[,,'Clinic']+tau[,,'Home'])
#     means[, ,'M'] <- (tau_combined*MDelta[,'M'] + norm_params['m_M']/norm_params['sigma_M']^2)/
#       (tau_combined + 1/norm_params['sigma_M']^2) + rnorm(num_indiv*num_samples)/sqrt(tau_combined + 1/norm_params['sigma_M']^2)
#     means[, ,'Delta'] <- (tau_combined*MDelta[,'Delta'] + norm_params['m_Delta']/norm_params['sigma_Delta']^2)/
#       (tau_combined + 1/norm_params['sigma_Delta']^2) + rnorm(num_indiv*num_samples)/sqrt(tau_combined + 1/norm_params['sigma_Delta']^2)
#     # For convenience, store the means for Home and clinic
#     means[, , 'Home'] <- means[, , 'M'] + means[, , 'Delta'] 
#     means[, , 'Clinic'] <- means[, , 'M'] - means[, , 'Delta']
#     # Now compute the density at this value of the individual parameters for the true distribution
#     logpi <- log(dgamma(tau[,,'Home'],shape = gamma_params['alpha','Home'], rate = gamma_params['beta','Home']))+
#       log(dgamma(tau[,,'Clinic'],shape = gamma_params['alpha','Clinic'], rate = gamma_params['beta','Clinic']))+
#       log(dnorm(means[,,'M'],mean = norm_params['m_M'],sd = norm_params['sigma_M']))+
#       log(dnorm(means[,,'Delta'],mean = norm_params['m_Delta'],sd = norm_params['sigma_Delta']))
#     # ... and for the distribution we really sampled from
#     logq <- log(dgamma(tau[,,'Home'],shape= gamma_params['alpha','Home']+(k-1)/2, 
#                            rate = gamma_params['beta','Home']+S2k[,'Home']))+
#       log(dgamma(tau[,,'Clinic'],shape= gamma_params['alpha','Clinic']+(k-1)/2, 
#                            rate = gamma_params['beta','Clinic']+S2k[,'Clinic']))+
#       log(dnorm(means[ , , 'M' ],
#                 mean = (tau_combined*MDelta[,'M'] + norm_params['m_M']/ norm_params['sigma_M']^2)/
#       (tau_combined + 1/norm_params['sigma_M']^2),
#       sd= 1/sqrt(tau_combined + 1/norm_params['sigma_M']^2)))+
#       log(dnorm(means[ , ,'Delta'],
#                 mean = (tau_combined*MDelta[,'Delta'] + norm_params['m_Delta']/norm_params['sigma_Delta']^2)/
#       (tau_combined + 1/norm_params['sigma_Delta']^2),
#       sd= 1/sqrt(tau_combined + 1/norm_params['sigma_Delta']^2)))
#     # Now multiply this importance ratio by the probability of the observations.
#     observ_probs <- array(0,c(num_indiv,num_samples,2),dimnames = list(NULL,NULL, BP_place_names))
#     for (BPplace in BP_place_names){
#             for (which_sim in seq_len(num_samples)){
#               observ_probs[ ,which_sim,BPplace] <-
#                 apply( pnorm( upper_lower[[BPplace]]$upper , means[ , which_sim, BPplace], 1/sqrt(tau[ , which_sim, BPplace]) ) -
#                          pnorm( upper_lower[[BPplace]]$lower , means[ , which_sim, BPplace], 1/sqrt(tau[ , which_sim, BPplace]) ), 1 , prod) # Probability for interval
#             }
#           }
#     log_likelihood_weighted <- exp(logpi - logq) * apply(log(observ_probs), c(1,2), sum )
#     list( parameters = list( tau = tau, means = means), densities= list(logpi = logpi, logq = logq), OP = observ_probs, LLW = log_likelihood_weighted )
# }

# Try out the diagnostics on the real data. First put the estimated parameters into the right form.
#ggplot(subset(all_s2k, s2k <20 & HC == 'Clinic'), aes(x=s2k, fill=RS)) + geom_histogram(alpha=0.2, position="identity")


```