R/Param_Fn.R

In aaronmberger/Geo_dGLMM_habitat: Package for conducting spatial estimation of CPUE index standardization (fisheries, ecology, etc).

Param_Fn <-
function( Version, DataList, RhoConfig=c("Beta1"=0,"Beta2"=0,"Epsilon1"=0,"Epsilon2"=0) ){
  # 
  if(Version%in%c("geo_index_v3e","geo_index_v3d","geo_index_v3c","geo_index_v3b","geo_index_v3a","geo_index_v2j","geo_index_v2i")){                                                                                                                                                                                                                                                                                                                                                                         
    Return = list("ln_H_input"=c(0,0), "beta1_t"=qlogis(0.01*0.99*tapply(ifelse(DataList$b_i>0,1,0),INDEX=factor(DataList$t_i,levels=1:DataList$n_t),FUN=mean)), "gamma1_j"=rep(0,DataList$n_j), "lambda1_k"=rep(0,DataList$n_k), "logetaE1"=0, "logetaO1"=0, "logkappa1"=0, "logsigmaV1"=log(1), "logsigmaVT1"=log(1), "nu1_v"=rep(0,DataList$n_v), "nu1_vt"=matrix(0,nrow=DataList$n_v,ncol=DataList$n_t), "Omegainput1_s"=rep(0,DataList$n_s), "Epsiloninput1_st"=matrix(0,nrow=DataList$n_s,ncol=DataList$n_t), "beta2_t"=log(tapply(ifelse(DataList$b_i>0,DataList$b_i/DataList$a_i,NA),INDEX=factor(DataList$t_i,levels=1:DataList$n_t),FUN=mean,na.rm=TRUE)), "gamma2_j"=rep(0,DataList$n_j), "lambda2_k"=rep(0,DataList$n_k), "logetaE2"=0, "logetaO2"=0, "logkappa2"=0, "logsigmaV2"=log(1), "logsigmaVT2"=log(1), "logSigmaM"=c(log(5),qlogis(0.8),log(2),log(5)), "nu2_v"=rep(0,DataList$n_v), "nu2_vt"=matrix(0,nrow=DataList$n_v,ncol=DataList$n_t), "Omegainput2_s"=rep(0,DataList$n_s), "Epsiloninput2_st"=matrix(0,nrow=DataList$n_s,ncol=DataList$n_t))
  }
  if(Version%in%c("geo_index_v3i","geo_index_v3h","geo_index_v3g","geo_index_v3f")){
    Return = list("ln_H_input"=c(0,0), "beta1_t"=qlogis(0.01*0.99*tapply(ifelse(DataList$b_i>0,1,0),INDEX=factor(DataList$t_i,levels=1:DataList$n_t),FUN=mean)), "gamma1_j"=rep(0,DataList$n_j), "lambda1_k"=rep(0,DataList$n_k), "logetaE1"=0, "logetaO1"=0, "logkappa1"=0, "logsigmaV1"=log(1), "logsigmaVT1"=log(1), "Beta_mean1"=0, "logsigmaB1"=log(1), "Beta_rho1"=0, "Epsilon_rho1"=0, "nu1_v"=rep(0,DataList$n_v), "nu1_vt"=matrix(0,nrow=DataList$n_v,ncol=DataList$n_t), "Omegainput1_s"=rep(0,DataList$n_s), "Epsiloninput1_st"=matrix(0,nrow=DataList$n_s,ncol=DataList$n_t), "beta2_t"=log(tapply(ifelse(DataList$b_i>0,DataList$b_i/DataList$a_i,NA),INDEX=factor(DataList$t_i,levels=1:DataList$n_t),FUN=mean,na.rm=TRUE)), "gamma2_j"=rep(0,DataList$n_j), "lambda2_k"=rep(0,DataList$n_k), "logetaE2"=0, "logetaO2"=0, "logkappa2"=0, "logsigmaV2"=log(1), "logsigmaVT2"=log(1), "Beta_mean2"=0, "logsigmaB2"=log(1), "Beta_rho2"=0, "Epsilon_rho2"=0, "logSigmaM"=c(log(5),qlogis(0.8),log(2),log(5)), "nu2_v"=rep(0,DataList$n_v), "nu2_vt"=matrix(0,nrow=DataList$n_v,ncol=DataList$n_t), "Omegainput2_s"=rep(0,DataList$n_s), "Epsiloninput2_st"=matrix(0,nrow=DataList$n_s,ncol=DataList$n_t))
  }
  if(Version%in%c("geo_index_v3k","geo_index_v3j")){
    Return = list("ln_H_input"=c(0,0), "hyperparameters_z"=rep(0,3), "beta1_t"=qlogis(0.01*0.99*tapply(ifelse(DataList$b_i>0,1,0),INDEX=factor(DataList$t_i,levels=1:DataList$n_t),FUN=mean)), "gamma1_j"=rep(0,DataList$n_j), "lambda1_k"=rep(0,DataList$n_k), "logetaE1"=0, "logetaO1"=0, "logkappa1"=0, "logsigmaV1"=log(1), "logsigmaVT1"=log(1), "Beta_mean1"=0, "logsigmaB1"=log(1), "Beta_rho1"=0, "Epsilon_rho1"=0, "nu1_v"=rep(0,DataList$n_v), "nu1_vt"=matrix(0,nrow=DataList$n_v,ncol=DataList$n_t), "Omegainput1_s"=rep(0,DataList$n_s), "Epsiloninput1_st"=matrix(0,nrow=DataList$n_s,ncol=DataList$n_t), "beta2_t"=log(tapply(ifelse(DataList$b_i>0,DataList$b_i/DataList$a_i,NA),INDEX=factor(DataList$t_i,levels=1:DataList$n_t),FUN=mean,na.rm=TRUE)), "gamma2_j"=rep(0,DataList$n_j), "lambda2_k"=rep(0,DataList$n_k), "logetaE2"=0, "logetaO2"=0, "logkappa2"=0, "logsigmaV2"=log(1), "logsigmaVT2"=log(1), "Beta_mean2"=0, "logsigmaB2"=log(1), "Beta_rho2"=0, "Epsilon_rho2"=0, "logSigmaM"=c(log(5),qlogis(0.8),log(2),log(5)), "nu2_v"=rep(0,DataList$n_v), "nu2_vt"=matrix(0,nrow=DataList$n_v,ncol=DataList$n_t), "Omegainput2_s"=rep(0,DataList$n_s), "Epsiloninput2_st"=matrix(0,nrow=DataList$n_s,ncol=DataList$n_t))
  }
  if(Version%in%c("geo_index_v3m","geo_index_v3l")){
    Return = list("ln_H_input"=c(0,0), "hyperparameters_z"=rep(0,3), "beta1_t"=qlogis(0.01*0.99*tapply(ifelse(DataList$b_i>0,1,0),INDEX=factor(DataList$t_i,levels=1:DataList$n_t),FUN=mean)), "gamma1_j"=rep(0,DataList$n_j), "gamma1_tp"=matrix(0,nrow=DataList$n_t,ncol=DataList$n_p), "lambda1_k"=rep(0,DataList$n_k), "logetaE1"=0, "logetaO1"=0, "logkappa1"=0, "logsigmaV1"=log(1), "logsigmaVT1"=log(1), "Beta_mean1"=0, "logsigmaB1"=log(1), "Beta_rho1"=0, "Epsilon_rho1"=0, "nu1_v"=rep(0,DataList$n_v), "nu1_vt"=matrix(0,nrow=DataList$n_v,ncol=DataList$n_t), "Omegainput1_s"=rep(0,DataList$n_s), "Epsiloninput1_st"=matrix(0,nrow=DataList$n_s,ncol=DataList$n_t), "beta2_t"=log(tapply(ifelse(DataList$b_i>0,DataList$b_i/DataList$a_i,NA),INDEX=factor(DataList$t_i,levels=1:DataList$n_t),FUN=mean,na.rm=TRUE)), "gamma2_j"=rep(0,DataList$n_j), "gamma2_tp"=matrix(0,nrow=DataList$n_t,ncol=DataList$n_p), "lambda2_k"=rep(0,DataList$n_k), "logetaE2"=0, "logetaO2"=0, "logkappa2"=0, "logsigmaV2"=log(1), "logsigmaVT2"=log(1), "Beta_mean2"=0, "logsigmaB2"=log(1), "Beta_rho2"=0, "Epsilon_rho2"=0, "logSigmaM"=c(log(5),qlogis(0.8),log(2),log(5)), "nu2_v"=rep(0,DataList$n_v), "nu2_vt"=matrix(0,nrow=DataList$n_v,ncol=DataList$n_t), "Omegainput2_s"=rep(0,DataList$n_s), "Epsiloninput2_st"=matrix(0,nrow=DataList$n_s,ncol=DataList$n_t))
  }
  # If either beta or epsilon is a random-walk process, fix starting value at 1
  if( "Beta_rho1"%in%names(Return) && RhoConfig[["Beta1"]]==2 ) Return[["Beta_rho1"]] = 1
  if( "Beta_rho2"%in%names(Return) && RhoConfig[["Beta2"]]==2 ) Return[["Beta_rho2"]] = 1
  if( "Epsilon_rho1"%in%names(Return) && RhoConfig[["Epsilon1"]]==2 ) Return[["Epsilon_rho1"]] = 1
  if( "Epsilon_rho2"%in%names(Return) && RhoConfig[["Epsilon2"]]==2 ) Return[["Epsilon_rho2"]] = 1
  # replace missing values function
  tmpfn = function( vec ){
    Return = ifelse( abs(vec)==Inf, NA, vec)
    return( ifelse(is.na(Return),mean(Return,na.rm=TRUE),Return) )
  }
  Return[["beta1_t"]] = tmpfn( Return[["beta1_t"]] )
  Return[["beta2_t"]] = tmpfn( Return[["beta2_t"]] )
  # Error messages
  if( any(sapply(Return, FUN=function(num){any(is.na(num))})) ) stop("Some parameter is NA")
  # Return tagged list
  return( Return )
}