`panel3pt8.fn` <-
function(){
# R wrapper for fitting occupancy model to the willow tit data using the likelihood specification given in Panel 3.8. This fits an occupancy model with covariates in both detection and occurrence probability using the zero-inflated binomial representation. Requires the willow tit data and utility functions
source("utilfns.Rd")
data<- read.table("wtmatrix.csv",header=TRUE,sep=",",na.strings=c("NA"))

elev<-data[,"elev"]
forest<-data[,"forest"]
elev<-scale(elev,center=TRUE)
forest<-scale(forest,center=TRUE)
ymat<-data[,c("y.1","y.2","y.3")]
M<-nrow(ymat)

length<-data[,"length"]
forest<-scale(data[,"forest"])
elev<-scale(data[,"elev"])
date<-as.matrix(data[,c("day.1","day.2","day.3")])

mdate<-mean(date,na.rm=TRUE)
sddate<-sqrt(var(date[1:length(date)],na.rm=TRUE))
date<-(date-mdate)/sddate

dur<-as.matrix(data[,c("dur.1","dur.2","dur.3")])
intensity<-dur/length
mint<-mean(intensity,na.rm=TRUE)
sdint<-sqrt(var(intensity[1:length(intensity)],na.rm=TRUE))
intensity<-(intensity-mint)/sdint
##length<- log(length)
length<- (length-mean(length))/sqrt(var(length))

lik<-function(parms,vars){
    tmp<-rep(0,9)
    names(tmp)<-c("pconst","psiconst","length","elev1","elev2","forest", "intensity","date1","date2")
    tmp[vars]<-parms
    ones<-rep(1,M)
    pmat<-expit(tmp[1]*ones + tmp[7]*intensity + tmp[8]*date + tmp[9]*(date^2))
    psi<-expit(tmp[2]*ones + tmp[3]*length + tmp[4]*elev + tmp[5]*(elev^2) + tmp[6]*forest)
    loglik<-rep(NA,M)
    for(i in 1:M){
       yvec<-ymat[i,]
       navec<-is.na(yvec)
       nd<-sum(yvec[!navec])
       pvec<-pmat[i,]
       cp<- (pvec^yvec)*((1-pvec)^(1-yvec))
       cp[navec]<-1
       loglik[i]<-log(prod(cp)*psi[i] + ifelse(nd==0,1,0)*(1-psi[i]))
   }
    sum(-1*loglik)
} 
nam<-c("pconst","psiconst","elev1","elev2","forest","length","intensity","date1","date2","aic")
out<-matrix(NA,nrow=10,ncol=length(nam))
dimnames(out)<-list(1:10,nam)

v<-c("pconst","psiconst","elev1","forest")
x<-nlm(lik,c(0,0,0,0),vars=v,hessian=TRUE)
out[1,v]<-x$estimate
out[1,"aic"]<-2*x$minimum + 2*length(v)

v<-c("pconst","psiconst","elev1","elev2","forest")
x<-nlm(lik,c(0,0,0,0,0),vars=v,hessian=TRUE)
out[2,v]<-x$estimate
out[2,"aic"]<-2*x$minimum + 2*length(v)


v<-c("pconst","psiconst","elev1","elev2","forest","length")
x<-nlm(lik,c(0,0,0,0,0,0),vars=v,hessian=TRUE)
out[3,v]<-x$estimate
out[3,"aic"]<-2*x$minimum + 2*length(v)

v<-c("pconst","psiconst","elev1","elev2","forest","length","intensity")
x<-nlm(lik,c(0,0,0,0,0,0,0),vars=v,hessian=TRUE)
out[4,v]<-x$estimate
out[4,"aic"]<-2*x$minimum + 2*length(v)

v<-c("pconst","psiconst","elev1","elev2","forest","length","intensity","date1")
x<-nlm(lik,c(0,0,0,0,0,0,0,0),vars=v,hessian=TRUE)
out[5,v]<-x$estimate
out[5,"aic"]<-2*x$minimum + 2*length(v)

v<-c("pconst","psiconst","elev1","elev2","forest","length","intensity","date1","date2")
x<-nlm(lik,c(0,0,0,0,0,0,0,0,0),vars=v,hessian=TRUE)
out[6,v]<-x$estimate
out[6,"aic"]<-2*x$minimum + 2*length(v)

v<-c("pconst","psiconst","intensity")
x<-nlm(lik,c(0,0,0),vars=v,hessian=TRUE)
out[7,v]<-x$estimate
out[7,"aic"]<-2*x$minimum + 2*length(v)

v<-c("pconst","psiconst","intensity","date1")
x<-nlm(lik,c(0,0,0,0),vars=v,hessian=TRUE)
out[8,v]<-x$estimate
out[8,"aic"]<-2*x$minimum + 2*length(v)

v<-c("pconst","psiconst","intensity","date1","date2")
x<-nlm(lik,c(0,0,0,0,0),vars=v,hessian=TRUE)
out[9,v]<-x$estimate
out[9,"aic"]<-2*x$minimum + 2*length(v)

v<-c("pconst","psiconst","length","intensity","date1","date2")
x<-nlm(lik,c(0,0,0,0,0,0),vars=v,hessian=TRUE)
out[10,v]<-x$estimate
out[10,"aic"]<-2*x$minimum + 2*length(v)



out 
 

}