R/dev_colextfp.r

In oliviergimenez/occuHMM: Fit occupancy data with hidden Markov models

#' Deviance of dynamic occupancy models with false positives
#'
#' For internal use. This function computes the -log(likelihood) of dynamic occupancy model with false positives.
#' @param
#' @keywords deviance
#' @export
#' @examples
#' dev_colextfp()

dev_colextfp <- function(b,data,eff,garb,nh,primary,secondary)
{
# b = vector of parameters on real scale
# data = site histories
# eff = nb of sites with that particular history
# garb = initial states
# nh = nb of sites
# primary/secondary occasions

#---------------------------------------------
# logit link
psi <- 1/(1+exp(-b[1]))
gamma <- 1/(1+exp(-b[2]))
epsilon <- 1/(1+exp(-b[3]))
p10 <- 1/(1+exp(-b[4]))
p11 <- 1/(1+exp(-b[5]))
delta <- 1/(1+exp(-b[6]))
#---------------------------------------------

#---------------------------------------------
# various quantities that will be useful later on
K <- length(primary)
J2 <- length(secondary)
J <- J2/K
N <- J * K
#---------------------------------------------

#---------------------------------------------
# initial states prob
PI <- c(1-psi,psi)

# transition prob
Aprimary <- matrix(c((1-gamma),gamma,epsilon,(1-epsilon)),nrow=2,ncol=2,byrow=T)
Asecondary <- matrix(c(1,0,0,1),nrow=2,ncol=2,byrow=T)
A <- array(NA,c(2,2,N))
for (j in primary) A[1:2,1:2,j] <- Aprimary
for (k in secondary[-primary]) A[1:2,1:2,k] <- Asecondary

# obs prob
B <- matrix(c(1-p10,1-p11,0,delta*p11,p10,(1-delta)*p11),nrow=3,ncol=2,byrow=T)
#---------------------------------------------

#---------------------------------------------
# calculate -log(lik)
l <- 0
for (i in 1:nh) # loop on sites
{
	oe <- garb[i] + 1 # first obs
	evennt <- data[,i] + 1 # non-det/det -> 1/2
	ALPHA <- PI*B[oe,]
	for (j in 2:N){ # loop on time
		ALPHA <- (ALPHA %*% A[1:2,1:2,j-1])*B[evennt[j],]
		}
      	l <- l + logprot(sum(ALPHA))*eff[i]
   }
l <- -l
l
#---------------------------------------------

}