"Copy" function to estimate surveillance sensitivity?

[Laura Cooper]

Hi INLA discussion group,

I have some surveillance data for an infectious disease in space and time. The observations consist of positive or negative results for two types of surveillance. I would like to estimate the sensitivity of each surveillance method (which we can call Y and Z) by estimating a latent space-time infection probability. I found some examples on state-space modelling in INLA using the "copy" function and I thought this might be an appropriate approach. However, I made a toy example with a single location to test and I am not getting the results I expected. I wonder if you can tell me what I am doing wrong?

I'll define my simple model first.

I have an infection probability i_t which is a first-order random walk:

# i is the hidden infection status
# i_t = a_t * i_t-1 + o_t
# o_t ~ N(0,O)
# Temporal coefficient a_t is constant and equal to one
# Therefore
# 0 = i_t - i_t-1 - o_t [this is a first-order random walk]

# I have some observations from surveillance Y_t and Z_t  with ny_t and nz_t number of tests which have a constant sensitivity beta_Y and beta_Z to detect i_t, (0<=beta_Y<=1 and 0<=beta_Z<=1).

# observation model
# Y_t ~ Binomial(p_t, ny_t)
# logit(p_t) = beta_Y * i_t
# Z_t ~ Binomial(q_t, nz_t)
# logit(q_t) = beta_Z * i_t

# Following this framework, I simulate some data:

library(INLA)

library(boot)

library(foreach)

i<-(sin(seq(0,3.14*4,by=0.05)-3.14*0.5)+1)/2
t<-1:length(i)

beta_Y<-0.8
beta_Z<-0.4

n_Y<-rep(20,times=length(t))
pos_Y<-unlist(foreach(x=t)%do%rbinom(1,n_Y[x],i[x]*beta_Y))

n_Z<-rep(20,times=length(t))
pos_Z<-unlist(foreach(x=t)%do%rbinom(1,n_Z[x],i[x]*beta_Z))

# I then set up my data with three response variables, with the first column representing Y_t, the second Z_t, and the third column 0 = i_t - i_t-1 - o_t.

# Number of time steps
n<-max(t)
# Because i_t depends on i_t-1
m<-n-1
# Response matrix (first column Y_t, second column Z_t, third column 0)
resp<-matrix(NA, n+n+m, 3)
resp[1:n, 1]<- pos_Y
resp[(n+1):(n+n), 2]<- pos_Z
resp[(n+n+1):(n+n+m),3]<- 0

# Three likelihoods
fam<-c("binomial","binomial","gaussian")

num_trials<-c(n_Y,n_Z,rep(NA,m))


i1 = c(1:n,1:n,2:n)             # indices for i_t
j = c(rep(NA,n),rep(NA,n),2:n -1)    #indices for i_t-1
w1 = c(rep(NA,n),rep(NA,n), rep(-1,m)) # weight for -1 * o_t-1
l = c(rep(NA,n),rep(NA,n),2:n -1) # indices for o_t
w2 <- c(rep(NA,n),rep(NA,n), rep(-1,m)) # weight for  -1 * o_(t-1)
i2 = c(1:n,rep(NA,n),rep(NA,m))             # indices for i_t method Y
i3 = c(rep(NA,n),1:n,rep(NA,m))           # indices for i_t method Z

# Then I fit my model and plot the results:

mod.inla <- inla(resp ~ 0 +
                   f(i1, model = "iid",hyper = list(prec = list(initial = -10, fixed = TRUE)))+
                   f(j, w1, copy = "i1") + f(l, w2, model = "iid")+
                   f(i2, copy = "i1",hyper = list(beta = list(fixed = FALSE))) +
                   f(i3, copy = "i1",hyper = list(beta = list(fixed = FALSE))) ,
                 data = list(resp=resp, i1=i1, i2=i2, i3=i3, j=j, l=l, w1=w1, w2=w2), family = fam,
                 control.family = list(list(), list(),
                                       list(hyper = list(prec = list(initial = 10, fixed = TRUE)))),
                 control.predictor = list(compute = TRUE),
                 Ntrials=num_trials
)

summary(mod.inla)

# Predicted response for Y compared to what I was expecting
matplot(mod.inla$summary.fitted.values[which(!is.na(resp[,1])),c(1,3,5)],col=2,type="l",lty=c(1,2,2),ylim=c(0,1))
lines(t,i*beta_Y)

# Predicted response for Z compared to what I was expecting

matplot(mod.inla$summary.fitted.values[which(!is.na(resp[,2])),c(1,3,5)],col=4,type="l",lty=c(1,2,2))
lines(t,i*beta_Z)

# Fitted i_t compared to true i_t

d<-mod.inla$summary.random$i1[,c(2,4,6)]
for(j in 1:3) d[,j]<-inv.logit(unlist(d[,j]))
matplot(d,col=3,type="l",lty=c(1,2,2),ylim=c(0,1))
lines(t,i)

# beta_Y and beta_Z (I think this is logit-linked?)
d<-mod.inla$summary.hyperpar[2:3,c(1,3,5)]
for(j in 1:3) d[,j]<-inv.logit(unlist(d[,j]))
d

# beta_Y compared to input
d<-inla.smarginal(mod.inla$marginals.hyperpar$`Beta for i2`)
d$x<-inv.logit(d$x)
plot(d,type="l",xlim=c(0,1))
abline(v=beta_Y,col=2)

# beta_Z compared to input
d<-inla.smarginal(mod.inla$marginals.hyperpar$`Beta for i3`)
d$x<-inv.logit(d$x)
plot(d,type="l",xlim=c(0,1))
abline(v=beta_Z,col=4)

Thanks in advance for your help and advice! Sorry if I have made some obvious mistakes - I am new to using INLA and have never used the "copy" functionality before.

Best,

Laura

[Helpdesk]

I think this one do what you try to do...

> --
> You received this message because you are subscribed to the Google
> Groups "R-inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to r-inla-discussion...@googlegroups.com.
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/83a1fc5e-d66f-460f-bf64-cffca71f2decn%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org