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)))