A proposal for fitting a model with multiple random effects

[JDChallenger]

#Simple GLMM. Based on mosquito death in a experimental hut trial
#Here random effects are crossed, rather than nested.
#But the overall method should apply more generally

#Can provide stan model comparison, if useful
library(ggplot2)
library(drjacoby)

#In this trial, we are evaluating two mosquito nets, one of which is an untreated (control) net
#There are four huts: each night a volunteer sleeps in each hut, under the net.
#The sleepers & the nets are moved around the different huts over the course of 32 nights

#Here's the trial structure. For 'net', 0=untreated control, 1=insecticide-treated net (ITN)
dfm <- data.frame('hut'=rep(c(1,2,3,4),32), 'sleeper' = c(rep(c(1,2,3,4),8), rep(c(4,1,2,3),8),
rep(c(3,4,1,2),8), rep(c(2,3,4,1),8)), 'net' = rep(c(0,0,1,1,1,1,0,0), 16) )
#How many mosquitoes enter each hut each night? Draw from a Poisson distribution
dfm$n <- rpois(128, lambda = 11)
dfm$dead <- NA

#To simulate how many mosquitoes die each night (binomial sampling), we include an effect for the ITN
#We also include a weaker effect, depending on the hut & sleeper
for(i in 1:128){
rho <- -1 + 0.1dfm$hut[i] - 0.2dfm$sleeper[i] + 1.5 * dfm$net[i]#prob on the logistic scale
#whack a bit o' noise on it
rho2 <- rnorm(1,mean = rho, sd = 0.1)
#convert to probability scale
prb <- exp(rho2)/(1+exp(rho2))
dfm$dead[i] <- rbinom(1, size = dfm$n[i], prob = prb)
}

#Now prepare what we need for drjacoby. There are 4 random effects for hut and 4 for sleeper.
#In total, there are 16 combinations of these
M <- 16 # no. of combinations of the two random effects? (4 values each)

df_params <- define_params(name = "hut1", min = -Inf, max = Inf, init = 0, block = c(1,2,3,4,M+1),
name = "hut2", min = -Inf, max = Inf, init = 0, block = c(5,6,7,8,M+1),
name = "hut3", min = -Inf, max = Inf, init = 0, block = c(9,10,11,12,M+1),
name = "hut4", min = -Inf, max = Inf, init = 0, block = c(13,14,15,16,M+1),
name = "sleeper1", min = -Inf, max = Inf, init = 0, block = c(1,5,9,13,M+1),
name = "sleeper2", min = -Inf, max = Inf, init = 0, block = c(2,6,10,14,M+1),
name = "sleeper3", min = -Inf, max = Inf, init = 0, block = c(3,7,11,15,M+1),
name = "sleeper4", min = -Inf, max = Inf, init = 0, block = c(4,8,12,16,M+1),
name = "sigma_hut", min = 0, max = Inf, init = 1, block = M+1,
name = "sigma_sleeper", min = 0, max = Inf, init = 1, block = M+1,
name = "a1", min = -Inf, max = Inf, init = 0, block = 1:M,
name = "a0", min = -Inf, max = Inf, init = 0, block = 1:M)

#Now we have to work out which data points need to go with each block
#empty list for the data
lst <- list()
#empty dataframe
dfa <- data.frame('blck'=rep(0,16), 'blckH'=rep(0,16), 'blckS'=rep(0,16))
count <- 1
for(j in 1:4){
for(k in 1:4){
aux <- dfm[dfm$hut == j & dfm$sleeper == k,]
#preare data for a block, store as a list
lst[[as.character(count)]] <- list(total = aux$n, dead = aux$dead, net = aux$net)
dfa$blck[count] <- count
dfa$blckH[count] <- j
dfa$blckS[count] <- k
count <- count + 1
}
}
#Look at dfa. For each block, which random effect do we need for hut (H) and sleeper (S)?
head(dfa)

#define log-likelihood function
r_loglike <- function(params, data, misc) {

#unpack parameters
#Patient-specific first
hut <- rep(0,4)
sleeper <- rep(0,4)
for (i in 1:4) {
hut[i] <- params[i]
sleeper[i] <- params[4 + i]
}
a1 <- params["a1"]
a0 <- params["a0"]
sigma_hut <- params["sigma_hut"]
sigma_sleeper <- params["sigma_sleeper"]

#get current update block
block = misc$block

#distinct method for first M blocks, vs (M+1)
ret = 0.0;
if (block == (M+1)) { # // likelihood for the global parameters
for(i in 1:4){
ret <- ret + dnorm(hut[i], mean=0, sd=sigma_hut, log=T) +
dnorm(sleeper[i], mean=0, sd=sigma_sleeper, log=T)
}
}else{
tot <- data[[block]]$total
ded <- data[[block]]$dead
nt <- data[[block]]$net
nnn <- length(tot)
for(j in 1:nnn){
#use dfa to index the random effects?
lprob <- a0 + a1*nt[j] + hut[dfa$blckH[block]] + sleeper[dfa$blckS[block]] # on logistic scale
#convert to probability scale
prb <- exp(lprob)/(1+exp(lprob))
ret <- ret + dbinom(ded[j], tot[j], prb, log = T)
}
}
return(ret)
}

#define log-prior function
r_logprior <- function(params, misc) {
sigma_hut = params["sigma_hut"]
sigma_sleeper = params["sigma_sleeper"]
a0 = params["a0"]
a1 = params["a1"]
l1 <- dnorm(x=a0, mean = 0, sd = 2.5, log = T)
l2 <- dnorm(x=a1, mean = 0, sd = 1, log = T)
l3 <- dexp(x = sigma_sleeper, rate = 1, log = T)
l4 <- dexp(x = sigma_hut, rate = 1, log = T)
ret <- l1 + l2 + l3 + l4
return(ret)
}

date()
mcmc <- run_mcmc(data = lst,
df_params = df_params,
loglike = r_loglike,
logprior = r_logprior,
burnin = 2e3,
samples = 6e3,
chains = 3,
#cluster = cl,
silent = F)
date()
plot_par(mcmc, show = c('a0','a1','sigma_hut','sigma_sleeper'), phase = "sampling")
plot_credible(mcmc, show = c('a0','a1','sigma_hut','sigma_sleeper'))
#random effect for sleeper should be (more-or-less) descending (see line 19)
plot_credible(mcmc, show = paste0('sleeper',seq(1,4)))
#random effect for hut should be (more-or-less) asscending (see line 19)
plot_credible(mcmc, show = paste0('hut',seq(1,4)))