`panel4pt3.fn` <-
function(){
# This is the R construction of the likelihood for the Royle-Nichols
# model with covariates on both detection and mean abundance
# (Panel 4.3). It uses the Swiss bird data (willow tit)
# The function is SLOW SLOW SLOW -- you can write a more efficient one
#
source("utilfns.Rd")

# read willow tit data
data<- read.table("wtmatrix.csv",header=TRUE,sep=",",na.strings=c("NA"))

# extract covariates and standardize
ymat<-as.matrix(data[,c("y.1","y.2","y.3")])
M<-nrow(ymat)
elev<-data[,"elev"]
forest<-data[,"forest"]
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
date2<-date*date
elev2<-elev*elev
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<- (length-mean(length))/sqrt(var(length))


Nmax<-100

# define the likelihood
RNlikcov<-function(parms,vars=c("a0","b0") ){
  ones<-rep(1,M)
  tmp<-c(0,0,0,0,0,0,0,0,0)
  nam<-c("b0","a0","intensity","date1","date2","elev1","elev2","forest","length")
  names(tmp)<-nam
  tmp[vars]<-parms
  b0<-tmp[1];a0<-tmp[2];a1<-tmp[3];a2<-tmp[4];  a3<-tmp[5];
  b1<-tmp[6];b2<-tmp[7];b3<-tmp[8];b4<-tmp[9]

  rvec<-expit(a0*ones + a1*intensity + a2*date + a3*date2)
  lamvec<-exp(b0*ones + b1*elev + b2*elev2 + b3*forest + b4*length)
  lik<-rep(NA,M)

  for(i in 1:M){
    gN<-dpois(0:Nmax,lamvec[i])
    gN<-gN/sum(gN)
    dvec<-ymat[i,]
    naflag<-is.na(dvec)
    PMAT<-  1-outer( (1-rvec[i,]),0:Nmax,"^")
    LIK<- t((PMAT^dvec)*((1-PMAT)^(1-dvec)))
    LIK[,naflag]<-1
    likvec<-apply(LIK,1,prod)
    lik[i]<-sum(likvec*gN)
  }
 -1*sum(log(lik))
}

vars<-out<-matrix(NA,nrow=8,ncol=10)

outnam<-c("b0","intensity","date1","date2","a0","elev1","elev2","forest","length","AIC")
dimnames(vars)<-dimnames(out)<-list(1:8,outnam)

cat("Null model",fill=TRUE)
v<-c("a0","b0")
print(a<-nlm(RNlikcov,c(.10,.10),vars=v,hessian=TRUE))
out[1,v]<-a$estimate
out[1,"AIC"]<-a$minimum*2 + 2*length(v)
vars[1,v]<-sqrt(diag(solve(a$hessian)))

cat("Null model + intensity",fill=TRUE)
v<-c("a0","b0","intensity")
print(a<-nlm(RNlikcov,c(.10,.10,.10),vars=v,hessian=TRUE))
out[2,v]<-a$estimate
out[2,"AIC"]<-a$minimum*2 + 2*length(v)
vars[2,v]<-sqrt(diag(solve(a$hessian)))

cat("Null model + intenstiy + DATE1",fill=TRUE)
v<-c("a0","b0","intensity","date1")
print(a<-nlm(RNlikcov,c(.10,.10,0,0),vars=v,hessian=TRUE))
out[3,v]<-a$estimate
out[3,"AIC"]<-a$minimum*2 + 2*length(v)
vars[3,v]<-sqrt(diag(solve(a$hessian)))

cat("Null model + intense + DATE1 + DATE2",fill=TRUE)
v<-c("a0","b0","intensity","date1","date2")
print(a<-nlm(RNlikcov,c(.10,.10,0,0,0),vars=v,hessian=TRUE))
out[4,v]<-a$estimate
out[4,"AIC"]<-a$minimum*2 + 2*length(v)
vars[4,v]<-sqrt(diag(solve(a$hessian)))

cat("Null model + intensity + DATE1 + DATE2 + LENGTH",fill=TRUE)
v<-c("a0","b0","intensity","date1","date2","length")
print(a<-nlm(RNlikcov,c(.10,.10,0,0,0,0),vars=v,hessian=TRUE))
out[5,v]<-a$estimate
out[5,"AIC"]<-a$minimum*2 + 2*length(v)
vars[5,v]<-sqrt(diag(solve(a$hessian)))


cat("Null model + intensity + DATE1 + DATE2 + LENGTH + FOREST",fill=TRUE)
v<-c("a0","b0","intensity","date1","date2","length","forest")
print(a<-nlm(RNlikcov,c(.10,.10,0,0,0,0,0),vars=v,hessian=TRUE))
out[6,v]<-a$estimate
out[6,"AIC"]<-a$minimum*2 + 2*length(v)
vars[6,v]<-sqrt(diag(solve(a$hessian)))

cat("Null model + intensity + DATE1 + DATE2 + LENGTH + FOREST + ELEV",fill=TRUE)
v<-c("a0","b0","intensity","date1","date2","length","forest","elev1")
print(a<-nlm(RNlikcov,c(.10,.10,0,0,0,0,0,0),vars=v,hessian=TRUE))
out[7,v]<-a$estimate
out[7,"AIC"]<-a$minimum*2 + 2*length(v)
vars[7,v]<-sqrt(diag(solve(a$hessian)))

cat("Null model + intensity + DATE1 + DATE2 + LENGTH + FOREST + ELEV + ELEV2",fill=TRUE)
v<-c("a0","b0","intensity","date1","date2","length","forest","elev1","elev2")
print(a<-nlm(RNlikcov,c(.10,.10,0,0,0,0,0,0,0),vars=v,hessian=TRUE))
out[8,v]<-a$estimate
out[8,"AIC"]<-a$minimum*2 + 2*length(v)
vars[8,v]<-sqrt(diag(solve(a$hessian)))


list(out=out,vars=vars)
 

}