Running a zero-altered Poisson (ZAP) model n R-inla

[valentina ndolo]

Hi,

I am trying to run a ZAP model in inla and I keep getting this error:

Error in inla(fZA.mesh1, family = c("zeroinflatedpoisson0", "binomial"),  : 

  Number of families 2 does not match number of response variables 1492

Here is my code below:

# 19.7 ZAP model with spatial correlation

# 19.7.1 Running the model in R-INLA

UG4$Species01  <- ifelse(UG4$total_deaths==0, 0, 1) # Presence/absence part

UG4$SpeciesPos <- ifelse(UG4$total_deaths > 0, UG4$total_deaths, NA) # count part

# Spatial index for count part

w1.pos.index <- inla.spde.make.index(name = 'wpos', n.spde  = spde1$n.spde)

# Spatial index for presence/absence part

w1.01.index <- inla.spde.make.index(name = 'w01', n.spde  = spde1$n.spde)

Xm <- model.matrix(~bio6.std + bio18.std + elevation.std +

                     distance.water.std +s.pH.std + s.calcium.std + s.carbon.std + 

                     s.water.std + mean.ndvi.std,

                   data = UG4)

N <- nrow(UG4)

X <- data.frame(Intercept  = rep(1, N), 

                bio6.std             = Xm[, 2],

                bio18.std            = Xm[, 3], 

                elevation.std        = Xm[, 4],

                distance.water.std   = Xm[, 5],

                s.pH.std             = Xm[, 6],

                s.calcium.std        = Xm[, 7],

                s.carbon.std         = Xm[, 8],

                s.water.std          = Xm[, 9],

                mean.ndvi.std        = Xm[, 10])

X01 <- data.frame(Intercept.01         = rep(1, N),

                  bio6.std.01             = Xm[, 2],

                  bio18.std.01            = Xm[, 3], 

                  elevation.std.01        = Xm[, 4],

                  distance.water.std.01   = Xm[, 5],

                  s.pH.std.01             = Xm[, 6],

                  s.calcium.std.01        = Xm[, 7],

                  s.carbon.std.01         = Xm[, 8],

                  s.water.std.01          = Xm[, 9],

                  mean.ndvi.std.01        = Xm[, 10])

# And this is the stack for the ZAP model

StackPos.mesh1 <- inla.stack(

  tag  = "FitPos",

  data = list(AllY = cbind(UG4$SpeciesPos, NA)),  

  A    = list(1, A1),                      

  effects = list(            

    Xpos = as.data.frame(X),

    wpos = w1.pos.index))

Stack01.mesh1 <- inla.stack(

  tag  = "Fit01",

  data = list(AllY = cbind(NA, UG4$Species01)),  

  A    = list(1, A1),                      

  effects = list(      

    X01 = as.data.frame(X01),

    w01 = w1.01.index))

Stack.ZA.mesh1 <- inla.stack(StackPos.mesh1, Stack01.mesh1)

# Specify the model formula

fZA.mesh1  <- AllY ~ -1 + Intercept + bio6.std + bio18.std + elevation.std +

  distance.water.std +s.pH.std + s.calcium.std + s.carbon.std + 

  s.water.std + mean.ndvi.std + 

  f(wpos, model = spde1) +

  

  Intercept.01 + bio6.std.01 + bio18.std.01 + elevation.std.01 +

  distance.water.std.01 +s.pH.std.01 + s.calcium.std.01 + s.carbon.std.01 + 

  s.water.std.01 + mean.ndvi.std.01 + 

  f(w01, model = spde1)

# Run model 

HyperZap <- list(theta = list(initial = -10, fixed = TRUE)) # For zero-truncated

ZAP.mesh1 <- inla(fZA.mesh1,

                  family = c("zeroinflatedpoisson0", "binomial"),  

                  control.family = list(list(hyper = HyperZap),

                                        list()),

                  data = inla.stack.data(Stack.ZA.mesh1),

                  control.compute = list(dic = TRUE, waic = TRUE),

                  control.predictor = list(

                    link = 1,

                    A = inla.stack.A(Stack.ZA.mesh1)))

Many thanks in advance.

Best regards,

Valentina.

[Elias T. Krainski]

Hi,

what is the output of

str(cbind(UG4$SpeciesPos, NA))

str(cbind(NA, UG4$Species01))

best regards,

Elias

--
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/9e5e0a9c-708f-4082-b75b-863cff4d3594n%40googlegroups.com.

[valentina ndolo]

Hi Elias,

At first I ran it and got 'NULL" results for both, but then I re-ran the entire code and got the following results:

> str(cbind(UG4$SpeciesPos, NA))

 int [1:746, 1:2] 1 8 2 1 3 7 2 1 1 1 ...

> str(cbind(NA, UG4$Species01))

 num [1:746, 1:2] NA NA NA NA NA NA NA NA NA NA ...

It seems like one line of code was skipped before. It now works perfectly fine.

Many thanks for the clue.

Kind regards,

Valentina.