R/lbm__habitat.R

In AtlanticR/lbm: Lattice-based models

lbm__habitat = function( p, dat, pa ) {
   #\\ this is the core engine of lbm .. localised space-time habiat modelling
 
  if (0) {
    if (!exists("nsims", p)) p$nsims = 5000
    if (!exists("habitat.threshold.quantile", p)) p$habitat.threshold.quantile = 0.05 # quantile at which to consider zero-valued abundance
  }

  sdTotal=sd(dat[,p$variables$Y], na.rm=T)

  if ( exists("lbm_local_model_distanceweighted", p) ) {
    if (p$lbm_local_model_distanceweighted) {
      Hmodel = try( gam( p$lbm_local_modelformula, data=dat, family=binomial(), weights=weights, optimizer=c("outer","optim")  ) )
    } else {
      Hmodel = try( gam( p$lbm_local_modelformula, data=dat, family=binomial(), optimizer=c("outer","optim")  ) )
    }
  } else {
      Hmodel = try( gam( p$lbm_local_modelformula, data=dat, family=binomial() ) )
  } 
  if ( "try-error" %in% class(Hmodel) ) return( NULL )

  dat$P = try( predict( Hmodel, newdata=dat, type="response", se.fit=FALSE ) ) 
  dat$Yhat = dat$P * dat$A

  rsq = cor( dat$Yhat, dat[,p$variables$Y], use="pairwise.complete.obs" )^2
  if (rsq < p$lbm_rsquared_threshold ) return(NULL)

  Hmodel.coef = mvtnorm::rmvnorm(p$nsims, coef(Hmodel), Hmodel$Vp, method="chol")
  rm( Hmodel); gc()
  Hsim =  predict(Hmodel, newdata=pa, type="lpmatrix") %*% t(Hmodel.coef) 
  rm( Hmodel.coef); gc()
  oops = which( is.na(Hsim) )
  if (length(oops) > 0)  Hsim[oops ] = 0  # assume to be zero
  
  pa$logitmean = apply( Hsim, 1, mean, na.rm=T )
  pa$logitsd = apply( Hsim, 1, sd, na.rm=T )

  Amodel.coef = mvtnorm::rmvnorm(p$nsims, coef(Amodel), Amodel$Vp, method="chol")
  rm(Amodel); gc()
  Asim =  predict(Amodel, newdata=pa, type="lpmatrix") %*% t(Amodel.coef) 
  rm (Amodel.coef); gc()
  oops = which( is.na(Asim) )
  if (length(oops) > 0)  Asim[oops ] = 0  # assume to be zero
  
  # Do not extrapolate: trim to XX% quantiles to be a little more conservative
  oopu =  which( Asim > p$qs[2] )
  if (length(oopu) > 0)  Asim[ oopu ] = p$qs[2]

  oopl =  which( Asim < p$qs[1]  )
  if (length(oopl) > 0)  Asim[ oopl ] = 0  # below detection limits

  Asim = Asim * Hsim  # Asim now becomes weighted by Pr of habitat

  pa$mean = as.vector( apply( Asim, 1, mean, na.rm=T ) )
  pa$sd =  as.vector( apply( Asim, 1, sd, na.rm=T ) )

  # iHabitat = which( pa$logitmean > p$habitat.threshold.quantile & (pa$logitmean - 2 * pa$logitsd) > 0 )

  ss = summary(hmod)
  lbm_stats = list( sdTotal=sdTotal, rsquared=rsq, ndata=nrow(dat) ) # must be same order as p$statsvars

  return( list( predictions=pa, lbm_stats=lbm_stats ) )  

}