Nothing
#########
# Compute MLE index of abundance
#########
ComputeMleIndices = function(Data, Model, FileName, Folder=NA, Weights="StrataAreas", StrataTable, Run=TRUE){
# Make folder
if(is.na(Folder)) Folder = paste(getwd(),"/",sep="")
# Attach stuff -- listed by search()
attach(Model$BUGSoutput$sims.list)
#attach(Data)
modelStructure = Model$modelStructure
Dist = Model$likelihood
# Estimate marginal means
Mat = matrix(NA,nrow=length(levels(year)),ncol=length(levels(strata)))
DataNew = data.frame(year=levels(year)[as.vector(row(Mat))], strata=levels(strata)[as.vector(col(Mat))])
DataNew = data.frame(DataNew, logeffort=rep(log(1),nrow(DataNew)), vesselYear=rep(999,nrow(DataNew)), ones.vec=rep(log(1),nrow(DataNew))) # I need to include ones.vec as zero for some reason
# Only run if there's no random effects and the distribution is either Gamma and Lognormal
if( (modelStructure$VesselYear.zeroTows%in%c("zero","fixed")) & (modelStructure$VesselYear.positiveTows%in%c("zero","fixed")) & (modelStructure$StrataYear.zeroTows%in%c("zero","fixed")) & (modelStructure$StrataYear.positiveTows%in%c("zero","fixed")) & (Dist=="gamma" | Dist=="lognormal") & Run==TRUE){
# Default formulae
FormulaPres = " ~ 0 + factor(year)"
if(nlevels(strata)>1) FormulaPres = paste(FormulaPres," + factor(strata)",sep="")
FormulaPos = " ~ 0 + factor(year)"
if(nlevels(strata)>1) FormulaPos = paste(FormulaPos," + factor(strata)",sep="")
# Modified formulae
if(modelStructure$StrataYear.zeroTows=="fixed" & nlevels(strata)>1){FormulaPres = paste(FormulaPres," + factor(strata):factor(year)",sep="")}
#if(modelStructure$VesselYear.zeroTows=="fixed"){FormulaPres = paste(FormulaPres," + factor(vesselYear)",sep="")}
#if(modelStructure$Catchability.zeroTows=="linear"){FormulaPres = paste(FormulaPres," + logeffort",sep="")}
# if(modelStructure$Catchability.zeroTows=="quadratic"){FormulaPres = paste(FormulaPres," + logeffort + logeffort2",sep="")}
if(modelStructure$StrataYear.positiveTows=="fixed" & nlevels(strata)>1){FormulaPos = paste(FormulaPos," + factor(strata):factor(year)",sep="")}
#if(modelStructure$VesselYear.positiveTows=="fixed"){FormulaPos = paste(FormulaPos," + factor(vesselYear)",sep="")}
#if(modelStructure$Catchability.positiveTows=="linear"){FormulaPos = paste(FormulaPos," + logeffort",sep="")}
# if(modelStructure$Catchability.positiveTows=="quadratic"){FormulaPos = paste(FormulaPos," + logeffort + logeffort2",sep="")}
#### Presence/absence
GlmPres <- glm(as.formula(paste("ifelse(Data[,'HAUL_WT_KG']>0,1,0)",FormulaPres)), family=binomial, control=glm.control(epsilon=1e-8, maxit=1000, trace=FALSE))
#### Positive catches
#OffsetPos = list(ones.vec,logeffort)[[ifelse(modelStructure$Catchability.positiveTows=="one",2,1)]]
if(Dist=="gamma") GlmPos <- glm(as.formula(paste("HAUL_WT_KG",FormulaPos)), family=Gamma(link="log"), offset=logeffort, subset=which(Data[,'HAUL_WT_KG']>0), control=glm.control(epsilon=1e-8, maxit=1000, trace=FALSE))
if(Dist=="lognormal") GlmPos <- lm(as.formula(paste("log(HAUL_WT_KG)",FormulaPos)), offset=logeffort, subset=which(Data[,'HAUL_WT_KG']>0))
# Calculate strata areas
Area = matrix(NA,nrow=length(unique(StrataTable[,'year'])),ncol=length(unique(StrataTable[,'strata'])))
for(YearI in 1:nrow(Area)){
for(StratI in 1:ncol(Area)){
AreaI = which(StrataTable[,'strataYear']==paste(levels(strata)[StratI],":",levels(year)[YearI],sep=""))
if(Weights=="StrataAreas") Area[YearI,StratI] = StrataTable[AreaI,'Area_Hectares']
if(Weights=="Equal") Area[YearI,StratI] = 1
}}
# Predict indices
IndexPres = array(predict(GlmPres, newdata=DataNew, type="response"), dim=dim(Mat))
IndexPos = array(predict(GlmPos, newdata=DataNew, type="response"), dim=dim(Mat))
if(Dist=="lognormal") IndexPos = exp(IndexPos + summary(GlmPos)$sigma^2/2) # Bias correction
Index = Area * IndexPres * IndexPos
}else{
Area = IndexPres = IndexPos = Index = matrix(NA,nrow=length(unique(StrataTable[,'year'])),ncol=length(unique(StrataTable[,'strata'])))
}
# EJ's code
#DglmData = data.frame(Catch=Data[,'HAUL_WT_KG']/Data[,'effort'], Year=Data[,'year'])
#Dglm = dglm(DglmData, dist="gamma", J=TRUE) # J=TRUE calculates jacknife estimates of CV
# Return results
Results1 = data.frame(year=levels(year)[as.vector(row(Mat))], strata=levels(strata)[as.vector(col(Mat))], Index=as.vector(Index), Pres=as.vector(IndexPres), Pos=as.vector(IndexPos))
Results2 = data.frame(year=levels(year), Index=rowSums(Index), Pres=rowSums(IndexPres), Pos=rowSums(IndexPos))
# Write and print output
write.csv(Results1,file=paste(Folder,"/",FileName,"ResultsByYearAndStrata_MLE.csv",sep=""))
write.csv(Results2,file=paste(Folder,"/",FileName,"ResultsByYear_MLE.csv",sep=""))
# Detach stuff -- listed by search()
#detach(Data)
detach(Model$BUGSoutput$sims.list)
return(list(Results1=Results1, Results2=Results2))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.