table8pt2.fn <-
function(){
### This function produces the results summarized in Table 8.2 
### (top panel).  The material in the bottom panel of Table 8.2
## can be produced using "panel4pt3.fn"

source("utilfns.Rd")
data<- read.table("wtmatrix.csv",header=TRUE,sep=",",na.strings=c("NA"))

elev<-data[,"elev"]
forest<-data[,"forest"]
elev<-as.numeric(scale(elev,center=TRUE))
forest<-as.numeric(scale(forest,center=TRUE))
ymat<-as.matrix(data[,c("c.1","c.2","c.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


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<- log(length)
length<- (length-mean(length))/sqrt(var(length))

Nmax<-100




likcov<-function(parms,vars=c("a0","b0") ){
ones<-rep(1,M)
tmp<-c(0,0,0,0,0,0,0,0,0)

nam<-c("a0","b0","intensity","date1","date2","elev1","elev2","forest","length")
names(tmp)<-nam
tmp[vars]<-parms

a0<-tmp[1]
b0<-tmp[2]
b1<-tmp[3]
b2<-tmp[4]
b3<-tmp[5]
a1<-tmp[6]
a2<-tmp[7]
a3<-tmp[8]
a4<-tmp[9]


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

LIK<-matrix(NA,nrow=Nmax+1,ncol=3)
for(i in 1:M){
    gN<-dpois(0:Nmax,lamvec[i])
    gN<-gN/sum(gN)
    dvec<-ymat[i,]
    naflag<-is.na(dvec)
    LIK[,1]<-dbinom(dvec[1],0:Nmax,rvec[i,1])
    LIK[,2]<-dbinom(dvec[2],0:Nmax,rvec[i,2])
    LIK[,3]<-dbinom(dvec[3],0:Nmax,rvec[i,3])
    LIK[,naflag]<-1
    likvec<-apply(LIK,1,prod)
    lik[i]<-sum(likvec*gN)
}

-1*sum(log(lik))
}






out<-matrix(NA,nrow=8,ncol=10)
vars<-out
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(likcov,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(likcov,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 + intensity + DATE1",fill=TRUE)
v<-c("a0","b0","intensity","date1")
print(a<-nlm(likcov,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 + intensity + DATE1 + DATE2",fill=TRUE)
v<-c("a0","b0","intensity","date1","date2")
print(a<-nlm(likcov,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(likcov,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(likcov,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(likcov,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(likcov,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)
 

}