Code for alli example

define.constr <- function(reference.i, reference.k, reference.j, alli,
                          ni = I, nk = K, nj = J)
{
    ## This function define the corner point constaint
    ##reference.i = reference category for  food
    ##reference.k = reference category for lake
    ##reference.j = reference category for size
   
    ## DEFINE CONSTRAINT   
    xx=numeric(ni+nk)
    kk=1
    for(k in 1:nk) {
        xx[kk]=idx.map(reference.i, k, ni)
        alli$ik[alli$ik==xx[kk]] = NA
        kk=kk+1
    }
    for(i in 1:ni) {
        xx[kk]=idx.map(i, reference.k, ni)
        alli$ik[alli$ik==xx[kk]] = NA
        kk=kk+1
    }
   
    xx=numeric(nj+nk)
    kk=1
    for(k in 1:nk) {
        xx[kk]=idx.map(reference.j, k, nj)
        alli$jk[alli$jk==xx[kk]] = NA
        kk=kk+1
    }
    for(j in 1:nj) {
        xx[kk]=idx.map(j, reference.k, nj)
        alli$jk[alli$jk==xx[kk]] = NA
        kk=kk+1
    }

    alli$food[alli$food==reference.k]=NA
    return(alli)
}

idx.map = function(i,j,ni)
{
    ## this function goes from the two index [i,j] to one index
    return (i + (j-1)*ni)
}

## Model with corner point constrains as in WINBUGS

alli = read.table("alli-dataset.dat", header=TRUE)

I = ni = 4 # number of food categories
J = nj = 2 # number of size categories
K = nk = 5 # number of lake categories


## get the index for beta_ik
nn =dim(alli)[1]
ik = numeric(nn)
ik.names=character(nn)
for(ii in 1:nn) {
    ik[ii] = idx.map(alli$lake[ii], alli$food[ii], ni)
    ik.names[ii] = paste("[",alli$lake[ii],",", alli$food[ii],"]",sep="")
}
alli$ik = ik
alli$ik.names = as.factor(ik.names)

## get the index for gamma_jk
jk = numeric(nn)
jk.names = character(nn)
for(ii in 1:nn) {
    jk[ii] = idx.map(alli$size[ii], alli$food[ii], nj)
    jk.names[ii] = paste("[",alli$size[ii],",", alli$food[ii],"]",sep="")

}
alli$jk = jk
alli$jk.names = as.factor(jk.names)

## get the index for lambda_ij (this is the extra parameter needed to
## go from a multinomial to the Poisson)
ij = numeric(nn)
for(ii in 1:nn) {
    ij[ii] = idx.map(alli$lake[ii], alli$size[ii], ni)
}
alli$ij = ij


## define the contraint here we use conrner point contraints as in the
## WinBUGS manual, the reference categories are as in WinBUGS

alli1 = define.constr(1,1,1,alli)
prior.prec = list(prec = list(initial = log(0.00001), fixed=TRUE))

## write the formula

formula = counts ~ f(ij, model="iid", hyper = prior.prec ) +
    f(food,  constr = FALSE, hyper = prior.prec) +
    f(ik, constr=FALSE, hyper = prior.prec) +
    f(jk, constr=FALSE, hyper = prior.prec) -1

## run the model here i have used the lincomb option to compute the
## results using the sum-to-zero constraint so that you can compare
## these resuts with those in alli_model1.R (is exactly the same thing
## done in winbugs)

beta12 = inla.make.lincomb(ik = c(-0.25,-0.25,-0.25))
names(beta12)="beta12"
beta13 = inla.make.lincomb(ik = c(rep(NA,3), -0.25, -0.25, -0.25))
names(beta13)="beta13"
beta14 = inla.make.lincomb(ik = c(rep(NA,6), -0.25, -0.25, -0.25))
names(beta14)="beta14"
beta15 = inla.make.lincomb(ik = c(rep(NA,9), -0.25, -0.25, -0.25))
names(beta15)="beta15"
all.lc = c(beta12, beta13, beta14, beta15)

mod1 = inla(formula, data = alli1, family = "poisson",
        control.predictor = list(compute=TRUE),
        lincomb = all.lc)




Comments