Code for model B of salamander data

###########################################
## Crossed effects - Salamander: Model B ##
###########################################

## load the INLA library
library (INLA)

## load the salamander data
load("salam.RData")

## organize data into a form suitable for logistic regression
dat0=data.frame("y"=c(salam$y),
    "fW"=as.integer(salam$x[,"W/R"]==1 | salam$x[,"W/W"]==1),
    "mW"=as.integer(salam$x[,"R/W"]==1 | salam$x[,"W/W"]==1),
    "WW"=as.integer(salam$x[,"W/W"]==1 ) )

## add salamander id (female, male)
id = t( apply(salam$z, 1, function(x) {
        tmp = which (x==1)
        tmp[2] = tmp[2] - 20
        tmp
    }) )

## index for the experiment group
dat0$group=as.factor(c(rep(c(rep(1,5),rep(2,10),rep(1,5)),6),
                       rep(c(rep(3,5),rep(4,10),rep(3,5)),6),
                       rep(c(rep(5,5),rep(6,10),rep(5,5)),6)))   
## index for the experiment
dat0$experiment=as.factor(rep(1:3, each=120))

## set the indices for the first two experiments modeled as iid2d,
## (The indices for the third experiment are set to NA)
fid_iid2_e1e2 <- c(id[,1],id[,1] + 20, rep(NA, 120))
mid_iid2_e1e2 <- c(id[,2],id[,2] + 20, rep(NA, 120))

## set the indices for third experiment
## (The indices for the first two experiments are set to NA)
fid_id_e3 <- c(rep(NA,240),id[,1])
mid_id_e3 <- c(rep(NA,240),id[,2])

## Indicator for fall
fall <- c(rep(0, 120), rep(1,240))

## generate the dataset
data <- data.frame(dat0, fid_iid2_e1e2, mid_iid2_e1e2, fid_id_e3, mid_id_e3, fall)


## specify the formula
formula <- y ~ fW+mW+WW+fall+
          f(fid_iid2_e1e2, model="iid2d", n=40, param=c(3,1.244,1.244,0)) +
          f(mid_iid2_e1e2, model="iid2d", n=40, param=c(3,1.244,1.244,0)) +
          f(fid_id_e3, model="iid",  param=c(1,0.622)) +
          f(mid_id_e3, model="iid",  param=c(1,0.622))
         
     
## Note: To get smooth posterior marginals for the hyperparameters you
## need to increase the number of maximum function evaluations:
## numint.maxfeval (default: 10000)
result <- inla(formula,data=data,family="binomial",Ntrials=rep(1,360),
  control.inla = list(numint.maxfeval = 200000), verbose=T)

summary(result)
plot(result)

Comments