occupancy models - Rota vs MacKenzie

gistfile1.txt

# we get same AIC when fitting Rota model with no interaction 
# vs fitting MacKenzie model to species separatately

library(unmarked)
library(mipfp)
?occuMulti

# simulate 2 species data
# https://github.com/oliviergimenez/2speciesoccupancy

set.seed(2022) 
S <- 2 # nb species 
N <- 500 # nb sites
J <- 5 # nb visits

psi01 <- 81/175
psi10 <- 36/175
psi11 <- 4/175
psi00 <- 1 - (psi01 + psi10 + psi11) # 54/175

psiS1 <- psi10 + psi11
psiS2 <- psi01 + psi11

or <- matrix(c(1, (psiS1*(1-psiS2))/(psiS2*(1-psiS1)), 
               (psiS2*(1-psiS1))/(psiS1*(1-psiS2)), 1), nrow = 2, ncol = 2, byrow = TRUE)
rownames(or) <- colnames(or) <- c("sp1", "sp2")
marg.probs <- c(psiS1, psiS2)
p.joint <- ObtainMultBinaryDist(odds = or, marg.probs = marg.probs)
z <- RMultBinary(n = N, mult.bin.dist = p.joint)$binary.sequences 
ps <- c(0.5,0.9)

y <- list()
for (i in 1:S){
  y[[i]] <- matrix(NA,N,J)
  for (j in 1:N){
    for (k in 1:J){
      y[[i]][j,k] <- rbinom(1,1,z[j,i]*ps[i])
    }
  }
}
names(y) <- c('sp1','sp2')

# fit Rota model with no interactions, get AIC

data <- unmarkedFrameOccuMulti(y=y)
occFormulas <- c('~1','~1','~0') 
detFormulas <- c('~1','~1')
fit <- occuMulti(detFormulas,occFormulas,data)
fit@AIC

# fit MacKenzie model to both species separately, and get sum of AIC

sp1 <- unmarkedFrameOccu(y=y[[1]])
fit1 <- occu(~1 ~1, sp1)
fit1
sp2 <- unmarkedFrameOccu(y=y[[2]])
fit2 <- occu(~1 ~1, sp2)
fit2

fit1@AIC + fit2@AIC