R/ComputeMleIndices.R

Defines functions ComputeMleIndices

Documented in ComputeMleIndices

#########
# 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))
}    

Try the nwfscDeltaGLM package in your browser

Any scripts or data that you put into this service are public.

nwfscDeltaGLM documentation built on May 2, 2019, 6:30 p.m.