From b8873a0c94aee88db13352ae68eaf2061871ab25 Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Wed, 12 Jul 2017 10:56:57 -0700 Subject: [PATCH] add V2.6.0... .... which allows independent factors for each category via input FieldConfig[X] = "IID" where X = {1,2,3,4} --- R/Data_Fn.R | 5 +- R/Make_Map.R | 18 +- R/Param_Fn.R | 10 +- inst/executables/VAST_v2_6_0.cpp | 983 +++++++++++++++++++++++++++++++ 4 files changed, 1004 insertions(+), 12 deletions(-) create mode 100644 inst/executables/VAST_v2_6_0.cpp diff --git a/R/Data_Fn.R b/R/Data_Fn.R index 14de6f4..4a64041 100644 --- a/R/Data_Fn.R +++ b/R/Data_Fn.R @@ -82,9 +82,10 @@ function( Version, FieldConfig, OverdispersionConfig=c("eta1"=0,"eta2"=0), ObsMo names(FieldConfig_input) = names(FieldConfig) g = function(vec) suppressWarnings(as.numeric(vec)) FieldConfig_input[] = ifelse( FieldConfig=="AR1", 0, FieldConfig_input) + FieldConfig_input[] = ifelse( FieldConfig=="IID", -2, FieldConfig_input) FieldConfig_input[] = ifelse( !is.na(g(FieldConfig)) & g(FieldConfig)>0 & g(FieldConfig)<=n_c, g(FieldConfig), FieldConfig_input) FieldConfig_input[] = ifelse( !is.na(g(FieldConfig)) & g(FieldConfig)==0, -1, FieldConfig_input) - if( any(is.na(FieldConfig_input)) ) stop( "'FieldConfig' must be: 0 (turn off overdispersion); 'AR1' (use AR1 structure); or 00 ){ List[[which(names(List)==list_names[1])]] = rnorm(n_f*n_c - n_f*(n_f-1)/2) List[[which(names(List)==list_names[2])]] = rarray(dim=as.vector(na.omit(c(n_i,n_f,n_t))), sd=sd) } + # AR1 process if( n_f==0 ){ List[[which(names(List)==list_names[1])]] = c(1,0.5) # Pointwise SD / Correlation List[[which(names(List)==list_names[2])]] = rarray(dim=as.vector(na.omit(c(n_i,n_c,n_t))), sd=sd) } + # Turned off if( n_f== -1 ){ List[[which(names(List)==list_names[1])]] = 1 # Turn off SD when zero factors, i.e., n_f = -1 List[[which(names(List)==list_names[2])]] = rarray(dim=as.vector(na.omit(c(n_i,abs(n_f),n_t))), sd=sd) } + # IID process + if( n_f== -2 ){ + List[[which(names(List)==list_names[1])]] = rep(1,n_c) # Pointwise SD / Correlation + List[[which(names(List)==list_names[2])]] = rarray(dim=as.vector(na.omit(c(n_i,n_c,n_t))), sd=sd) + } return( List ) } # @@ -78,7 +86,7 @@ function( Version, DataList, RhoConfig=c("Beta1"=0,"Beta2"=0,"Epsilon1"=0,"Epsil if(Version%in%c("VAST_v1_9_0","VAST_v1_8_0","VAST_v1_7_0","VAST_v1_6_0")){ Return = list("ln_H_input"=c(0,0), "beta1_ct"=NA, "gamma1_j"=rep(0,DataList$n_j), "gamma1_ctp"=array(0,dim=c(DataList$n_c,DataList$n_t,DataList$n_p)), "lambda1_k"=rep(0,DataList$n_k), "L1_z"=NA, "L_omega1_z"=NA, "L_epsilon1_z"=NA, "logkappa1"=log(0.9), "Beta_mean1"=0, "logsigmaB1"=log(1), "Beta_rho1"=0, "Epsilon_rho1"=0, "eta1_vf"=NA, "Omegainput1_sf"=NA, "Epsiloninput1_sft"=NA, "beta2_ct"=NA, "gamma2_j"=rep(0,DataList$n_j), "gamma2_ctp"=array(0,dim=c(DataList$n_c,DataList$n_t,DataList$n_p)), "lambda2_k"=rep(0,DataList$n_k), "L2_z"=NA, "L_omega2_z"=NA, "L_epsilon2_z"=NA, "logkappa2"=log(0.9), "Beta_mean2"=0, "logsigmaB2"=log(1), "Beta_rho2"=0, "Epsilon_rho2"=0, "logSigmaM"=rep(1,DataList$n_c)%o%c(log(5),log(2),log(5)), "eta2_vf"=NA, "Omegainput2_sf"=NA, "Epsiloninput2_sft"=NA ) } - if(Version%in%c("VAST_v2_5_0","VAST_v2_4_0","VAST_v2_3_0","VAST_v2_2_0","VAST_v2_1_0","VAST_v2_0_0")){ + if(Version%in%c("VAST_v2_6_0","VAST_v2_5_0","VAST_v2_4_0","VAST_v2_3_0","VAST_v2_2_0","VAST_v2_1_0","VAST_v2_0_0")){ Return = list("ln_H_input"=c(0,0), "beta1_ct"=NA, "gamma1_j"=rep(0,DataList$n_j), "gamma1_ctp"=array(0,dim=c(DataList$n_c,DataList$n_t,DataList$n_p)), "lambda1_k"=rep(0,DataList$n_k), "L1_z"=NA, "L_omega1_z"=NA, "L_epsilon1_z"=NA, "logkappa1"=log(0.9), "Beta_mean1"=0, "logsigmaB1"=log(1), "Beta_rho1"=0, "Epsilon_rho1"=0, "log_sigmaratio1_z"=rep(0,ncol(DataList$t_iz)), "eta1_vf"=NA, "Omegainput1_sf"=NA, "Epsiloninput1_sft"=NA, "beta2_ct"=NA, "gamma2_j"=rep(0,DataList$n_j), "gamma2_ctp"=array(0,dim=c(DataList$n_c,DataList$n_t,DataList$n_p)), "lambda2_k"=rep(0,DataList$n_k), "L2_z"=NA, "L_omega2_z"=NA, "L_epsilon2_z"=NA, "logkappa2"=log(0.9), "Beta_mean2"=0, "logsigmaB2"=log(1), "Beta_rho2"=0, "Epsilon_rho2"=0, "log_sigmaratio2_z"=rep(0,ncol(DataList$t_iz)), "logSigmaM"=rep(1,DataList$n_c)%o%c(log(5),log(2),log(5)), "eta2_vf"=NA, "Omegainput2_sf"=NA, "Epsiloninput2_sft"=NA ) } # Overdispersion diff --git a/inst/executables/VAST_v2_6_0.cpp b/inst/executables/VAST_v2_6_0.cpp new file mode 100644 index 0000000..0d5c6a6 --- /dev/null +++ b/inst/executables/VAST_v2_6_0.cpp @@ -0,0 +1,983 @@ +#include +#include + +// Function for detecting NAs +template +bool isNA(Type x){ + return R_IsNA(asDouble(x)); +} + +// Posfun +template +Type posfun(Type x, Type lowerlimit, Type &pen){ + pen += CppAD::CondExpLt(x,lowerlimit,Type(0.01)*pow(x-lowerlimit,2),Type(0)); + return CppAD::CondExpGe(x,lowerlimit,x,lowerlimit/(Type(2)-x/lowerlimit)); +} + +// Variance +template +Type var( array vec ){ + Type vec_mod = vec - (vec.sum()/vec.size()); + Type res = pow(vec_mod, 2).sum() / vec.size(); + return res; +} + +// dlnorm +template +Type dlnorm(Type x, Type meanlog, Type sdlog, int give_log=0){ + //return 1/(sqrt(2*M_PI)*sd)*exp(-.5*pow((x-mean)/sd,2)); + Type logres = dnorm( log(x), meanlog, sdlog, true) - log(x); + if(give_log) return logres; else return exp(logres); +} + +// Generate loadings matrix +template +matrix loadings_matrix( vector L_val, int n_rows, int n_cols ){ + matrix L_rc(n_rows, n_cols); + int Count = 0; + for(int r=0; r=c){ + L_rc(r,c) = L_val(Count); + Count++; + }else{ + L_rc(r,c) = 0.0; + } + }} + return L_rc; +} + +// IN: eta1_vf; L1_z +// OUT: jnll_comp; eta1_vc +template +matrix overdispersion_by_category_nll( int n_f, int n_v, int n_c, matrix eta_vf, vector L_z, Type &jnll_pointer){ + using namespace density; + matrix eta_vc(n_v, n_c); + vector Tmp_c; + // Turn off + if(n_f<0){ + eta_vc.setZero(); + } + // AR1 structure + if( n_f==0 ){ + for(int v=0; v0 ){ + // Assemble the loadings matrix + matrix L_cf = loadings_matrix( L_z, n_c, n_f ); + // Multiply out overdispersion + eta_vc = eta_vf * L_cf.transpose(); + // Probability of overdispersion + for(int v=0; v // +matrix convert_upper_cov_to_cor( matrix cov ){ + int nrow = cov.row(0).size(); + for( int i=0; i // +matrix gmrf_by_category_nll( int n_f, int method, int n_s, int n_c, Type logkappa, array gmrf_input_sf, array gmrf_mean_sf, vector L_z, density::GMRF_t gmrf_Q, Type &jnll_pointer){ + using namespace density; + matrix gmrf_sc(n_s, n_c); + Type logtau; + // IID + if(n_f == -2){ + for( int c=0; c0){ + if(method==0) logtau = log( 1 / (exp(logkappa) * sqrt(4*M_PI)) ); + if(method==1) logtau = log( 1 / sqrt(1-exp(logkappa*2)) ); + for( int f=0; f L_cf = loadings_matrix( L_z, n_c, n_f ); + gmrf_sc = (gmrf_input_sf.matrix() * L_cf.transpose()) / exp(logtau); + } + return gmrf_sc; +} + +// CMP distribution +template +Type dCMP(Type x, Type mu, Type nu, int give_log=0, int iter_max=30, int break_point=10){ + // Explicit + Type ln_S_1 = nu*mu - ((nu-1)/2)*log(mu) - ((nu-1)/2)*log(2*M_PI) - 0.5*log(nu); + // Recursive + vector S_i(iter_max); + S_i(0) = 1; + for(int i=1; i +Type dPoisGam( Type x, Type shape, Type scale, Type intensity, vector &diag_z, int maxsum=50, int minsum=1, int give_log=0 ){ + // Maximum integration constant to prevent numerical overflow, but capped at value for maxsum to prevent numerical underflow when subtracting by a higher limit than is seen in the sequence + Type max_log_wJ, z1, maxJ_bounded; + if( x==0 ){ + diag_z(0) = 1; + max_log_wJ = 0; + diag_z(1) = 0; + }else{ + z1 = log(intensity) + shape*log(x/scale) - shape*log(shape) + 1; + diag_z(0) = exp( (z1 - 1) / (1 + shape) ); + maxJ_bounded = CppAD::CondExpGe(diag_z(0), Type(maxsum), Type(maxsum), diag_z(0)); + max_log_wJ = maxJ_bounded*log(intensity) + (maxJ_bounded*shape)*log(x/scale) - lgamma(maxJ_bounded+1) - lgamma(maxJ_bounded*shape); + diag_z(1) = diag_z(0)*log(intensity) + (diag_z(0)*shape)*log(x/scale) - lgamma(diag_z(0)+1) - lgamma(diag_z(0)*shape); + } + // Integration constant + Type W = 0; + Type log_w_j, pos_penalty; + for( int j=minsum; j<=maxsum; j++ ){ + Type j2 = j; + //W += pow(intensity,j) * pow(x/scale, j2*shape) / exp(lgamma(j2+1)) / exp(lgamma(j2*shape)) / exp(max_log_w_j); + log_w_j = j2*log(intensity) + (j2*shape)*log(x/scale) - lgamma(j2+1) - lgamma(j2*shape); + //W += exp( posfun(log_w_j, Type(-30), pos_penalty) ); + W += exp( log_w_j - max_log_wJ ); + if(j==minsum) diag_z(2) = log_w_j; + if(j==maxsum) diag_z(3) = log_w_j; + } + // Loglikelihood calculation + Type loglike = 0; + if( x==0 ){ + loglike = -intensity; + }else{ + loglike = -x/scale - intensity - log(x) + log(W) + max_log_wJ; + } + // Return + if(give_log) return loglike; else return exp(loglike); +} + + +// Space time +template +Type objective_function::operator() () +{ + using namespace R_inla; + using namespace Eigen; + using namespace density; + + // Dimensions + DATA_INTEGER(n_i); // Number of observations (stacked across all years) + DATA_INTEGER(n_s); // Number of "strata" (i.e., vectices in SPDE mesh) + DATA_INTEGER(n_x); // Number of real "strata" (i.e., k-means locations) + DATA_INTEGER(n_t); // Number of time-indices + DATA_INTEGER(n_c); // Number of categories (e.g., length bins) + DATA_INTEGER(n_j); // Number of static covariates + DATA_INTEGER(n_p); // Number of dynamic covariates + DATA_INTEGER(n_k); // Number of catchability variables + DATA_INTEGER(n_v); // Number of tows/vessels (i.e., levels for the factor explaining overdispersion) + DATA_INTEGER(n_l); // Number of indices to post-process + DATA_INTEGER(n_m); // Number of range metrics to use (probably 2 for Eastings-Northings) + + // Config + DATA_IVECTOR( Options_vec ); + // Slot 0 -- Aniso: 0=No, 1=Yes + // Slot 1 -- DEPRECATED + // Slot 2 -- AR1 on betas (year intercepts) to deal with missing years: 0=No, 1=Yes + // Slot 3 -- DEPRECATED + // Slot 4 -- DEPRECATED + // Slot 5 -- Upper limit constant of integration calculation for infinite-series density functions (Conway-Maxwell-Poisson and Tweedie) + // Slot 6 -- Breakpoint in CMP density function + // Slot 7 -- Whether to use SPDE or 2D-AR1 hyper-distribution for spatial process: 0=SPDE; 1=2D-AR1 + DATA_IVECTOR(FieldConfig); // Input settings (vector, length 4) + DATA_IVECTOR(OverdispersionConfig); // Input settings (vector, length 2) + DATA_IVECTOR(ObsModel); // Observation model + // Slot 0: Distribution for positive catches + // Slot 1: Link function for encounter probabilities + DATA_IVECTOR(Options); // Reporting options + // Slot 0: Calculate SD for Index_xctl + // Slot 1: Calculate SD for log(Index_xctl) + // Slot 2: Calculate mean_Z_ctm (i.e., center-of-gravity) + // Slot 3: DEPRECATED + // Slot 4: Calculate mean_D_tl and effective_area_tl + // Slot 5: Calculate standard errors for Covariance and Correlation among categories using factor-analysis parameterization + // Slot 6: Calculate synchrony for different periods specified via yearbounds_zz + // Slot 7: Calculate coherence and variance for Epsilon1_sct and Epsilon2_sct + DATA_IMATRIX(yearbounds_zz); + // Two columns, and 1+ rows, specifying first and last t for each period used in calculating synchrony + + // Data vectors + DATA_VECTOR(b_i); // Response (biomass) for each observation + DATA_VECTOR(a_i); // Area swept for each observation (km^2) + DATA_IVECTOR(c_i); // Category for each observation + DATA_IVECTOR(s_i); // Station for each observation + DATA_IMATRIX(t_iz); // Time-indices (year, season, etc.) for each observation + DATA_IVECTOR(v_i); // tows/vessels for each observation (level of factor representing overdispersion) + DATA_VECTOR(PredTF_i); // vector indicating whether an observatino is predictive (1=used for model evaluation) or fitted (0=used for parameter estimation) + DATA_MATRIX(a_xl); // Area for each "real" stratum(km^2) in each stratum + DATA_MATRIX(X_xj); // Covariate design matrix (strata x covariate) + DATA_ARRAY(X_xtp); // Covariate design matrix (strata x covariate) + DATA_MATRIX(Q_ik); // Catchability matrix (observations x variable) + DATA_IMATRIX(t_yz); // Matrix for time-indices of calculating outputs (abundance index and "derived-quantity") + DATA_MATRIX(Z_xm); // Derived quantity matrix + + // SPDE objects + DATA_STRUCT(spde,spde_t); + + // Aniso objects + DATA_STRUCT(spde_aniso,spde_aniso_t); + + // Sparse matrices for precision matrix of 2D AR1 process + // Q = M0*(1+rho^2)^2 + M1*(1+rho^2)*(-rho) + M2*rho^2 + DATA_SPARSE_MATRIX(M0); + DATA_SPARSE_MATRIX(M1); + DATA_SPARSE_MATRIX(M2); + + // Parameters + PARAMETER_VECTOR(ln_H_input); // Anisotropy parameters + // -- presence/absence + PARAMETER_MATRIX(beta1_ct); // Year effect + PARAMETER_VECTOR(gamma1_j); // Static covariate effect + PARAMETER_ARRAY(gamma1_ctp); // Dynamic covariate effect + PARAMETER_VECTOR(lambda1_k); // Catchability coefficients + PARAMETER_VECTOR(L1_z); // Overdispersion parameters + PARAMETER_VECTOR(L_omega1_z); + PARAMETER_VECTOR(L_epsilon1_z); + PARAMETER(logkappa1); + PARAMETER(Beta_mean1); // mean-reversion for beta1_t + PARAMETER(logsigmaB1); // SD of beta1_t (default: not included in objective function) + PARAMETER(Beta_rho1); // AR1 for positive catch Epsilon component, Default=0 + PARAMETER(Epsilon_rho1); // AR1 for presence/absence Epsilon component, Default=0 + PARAMETER_VECTOR(log_sigmaratio1_z); // Ratio of variance for columns of t_iz + // -- Gaussian random fields + PARAMETER_MATRIX(eta1_vf); + PARAMETER_ARRAY(Omegainput1_sf); // Expectation + PARAMETER_ARRAY(Epsiloninput1_sft); // Annual variation + // -- positive catch rates + PARAMETER_MATRIX(beta2_ct); // Year effect + PARAMETER_VECTOR(gamma2_j); // Covariate effect + PARAMETER_ARRAY(gamma2_ctp); // Dynamic covariate effect + PARAMETER_VECTOR(lambda2_k); // Catchability coefficients + PARAMETER_VECTOR(L2_z); // Overdispersion parameters + PARAMETER_VECTOR(L_omega2_z); + PARAMETER_VECTOR(L_epsilon2_z); + PARAMETER(logkappa2); + PARAMETER(Beta_mean2); // mean-reversion for beta2_t + PARAMETER(logsigmaB2); // SD of beta2_t (default: not included in objective function) + PARAMETER(Beta_rho2); // AR1 for positive catch Epsilon component, Default=0 + PARAMETER(Epsilon_rho2); // AR1 for positive catch Epsilon component, Default=0 + PARAMETER_VECTOR(log_sigmaratio2_z); // Ratio of variance for columns of t_iz + PARAMETER_ARRAY(logSigmaM); // Slots: 0=mix1 CV, 1=prob-of-mix1, 2= + // -- Gaussian random fields + PARAMETER_MATRIX(eta2_vf); + PARAMETER_ARRAY(Omegainput2_sf); // Expectation + PARAMETER_ARRAY(Epsiloninput2_sft); // Annual variation + + // Indices -- i=Observation; j=Covariate; v=Vessel; t=Year; s=Stratum + int i,j,v,t,s,c,y,z; + + // Objective function + vector jnll_comp(13); + // Slot 0 -- spatial, encounter + // Slot 1 -- spatio-temporal, encounter + // Slot 2 -- spatial, positive catch + // Slot 3 -- spatio-temporal, positive catch + // Slot 4 -- tow/vessel overdispersion, encounter + // Slot 5 -- tow/vessel overdispersion, positive catch + // Slot 8 -- penalty on beta, encounter + // Slot 9 -- penalty on beta, positive catch + // Slot 10 -- likelihood of data, encounter + // Slot 11 -- likelihood of data, positive catch + jnll_comp.setZero(); + Type jnll = 0; + + // Derived parameters + Type Range_raw1, Range_raw2; + if( Options_vec(7)==0 ){ + Range_raw1 = sqrt(8) / exp( logkappa1 ); // Range = approx. distance @ 10% correlation + Range_raw2 = sqrt(8) / exp( logkappa2 ); // Range = approx. distance @ 10% correlation + } + if( Options_vec(7)==1 ){ + Range_raw1 = log(0.1) / logkappa1; // Range = approx. distance @ 10% correlation + Range_raw2 = log(0.1) / logkappa2; // Range = approx. distance @ 10% correlation + } + array SigmaM( n_c, 3 ); + SigmaM = exp( logSigmaM ); + + // Anisotropy elements + matrix H(2,2); + H(0,0) = exp(ln_H_input(0)); + H(1,0) = ln_H_input(1); + H(0,1) = ln_H_input(1); + H(1,1) = (1+ln_H_input(1)*ln_H_input(1)) / exp(ln_H_input(0)); + + //////////////////////// + // Calculate joint likelihood + //////////////////////// + + // Random field probability + Eigen::SparseMatrix Q1; + Eigen::SparseMatrix Q2; + GMRF_t gmrf_Q; + if( Options_vec(7)==0 & Options_vec(0)==0 ){ + Q1 = Q_spde(spde, exp(logkappa1)); + Q2 = Q_spde(spde, exp(logkappa2)); + } + if( Options_vec(7)==0 & Options_vec(0)==1 ){ + Q1 = Q_spde(spde_aniso, exp(logkappa1), H); + Q2 = Q_spde(spde_aniso, exp(logkappa2), H); + } + if( Options_vec(7)==1 ){ + Q1 = M0*pow(1+exp(logkappa1*2),2) + M1*(1+exp(logkappa1*2))*(-exp(logkappa1)) + M2*exp(logkappa1*2); + Q2 = M0*pow(1+exp(logkappa2*2),2) + M1*(1+exp(logkappa2*2))*(-exp(logkappa2)) + M2*exp(logkappa2*2); + } + // Probability of encounter + gmrf_Q = GMRF(Q1); + // Omega1 + int n_f; + n_f = Omegainput1_sf.cols(); + array Omegamean1_sf(n_s, n_f); + Omegamean1_sf.setZero(); + // Epsilon1 + n_f = Epsiloninput1_sft.col(0).cols(); + array Epsilonmean1_sf(n_s, n_f); + array Omega1_sc(n_s, n_c); + Omega1_sc = gmrf_by_category_nll(FieldConfig(0), Options_vec(7), n_s, n_c, logkappa1, Omegainput1_sf, Omegamean1_sf, L_omega1_z, gmrf_Q, jnll_comp(0)); + array Epsilon1_sct(n_s, n_c, n_t); + for(t=0; t=1){ + Epsilonmean1_sf = Epsilon_rho1 * Epsiloninput1_sft.col(t-1); + Epsilon1_sct.col(t) = gmrf_by_category_nll(FieldConfig(1), Options_vec(7), n_s, n_c, logkappa1, Epsiloninput1_sft.col(t), Epsilonmean1_sf, L_epsilon1_z, gmrf_Q, jnll_comp(1)); + } + } + // Positive catch rate + gmrf_Q = GMRF(Q2); + // Omega2 + n_f = Omegainput2_sf.cols(); + array Omegamean2_sf(n_s, n_f); + Omegamean2_sf.setZero(); + // Epsilon2 + n_f = Epsiloninput2_sft.col(0).cols(); + array Epsilonmean2_sf(n_s, n_f); + array Omega2_sc(n_s, n_c); + Omega2_sc = gmrf_by_category_nll(FieldConfig(2), Options_vec(7), n_s, n_c, logkappa2, Omegainput2_sf, Omegamean2_sf, L_omega2_z, gmrf_Q, jnll_comp(2)); + array Epsilon2_sct(n_s, n_c, n_t); + for(t=0; t=1){ + Epsilonmean2_sf = Epsilon_rho2 * Epsiloninput2_sft.col(t-1); + Epsilon2_sct.col(t) = gmrf_by_category_nll(FieldConfig(3), Options_vec(7), n_s, n_c, logkappa2, Epsiloninput2_sft.col(t), Epsilonmean2_sf, L_epsilon2_z, gmrf_Q, jnll_comp(3)); + } + } + + ////// Probability of correlated overdispersion among bins + // IN: eta1_vf; n_f; L1_z + // OUT: jnll_comp; eta1_vc + matrix eta1_vc(n_v, n_c); + eta1_vc = overdispersion_by_category_nll( OverdispersionConfig(0), n_v, n_c, eta1_vf, L1_z, jnll_comp(4) ); + matrix eta2_vc(n_v, n_c); + eta2_vc = overdispersion_by_category_nll( OverdispersionConfig(1), n_v, n_c, eta2_vf, L2_z, jnll_comp(5) ); + + // Possible structure on betas + if( Options_vec(2)!=0 ){ + for(c=0; c eta1_x = X_xj * gamma1_j.matrix(); + vector zeta1_i = Q_ik * lambda1_k.matrix(); + vector eta2_x = X_xj * gamma2_j.matrix(); + vector zeta2_i = Q_ik * lambda2_k.matrix(); + array eta1_xct(n_x, n_c, n_t); + array eta2_xct(n_x, n_c, n_t); + eta1_xct.setZero(); + eta2_xct.setZero(); + for(int x=0; x var_i(n_i); + Type var_y; + // Linear predictor (pre-link) for presence/absence component + vector P1_i(n_i); + // Response predictor (post-link) + // ObsModel = 0:4 or 11:12: probability ("phi") that data is greater than zero + // ObsModel = 5 (ZINB): phi = 1-ZeroInflation_prob -> Pr[D=0] = NB(0|mu,var)*phi + (1-phi) -> Pr[D>0] = phi - NB(0|mu,var)*phi + vector R1_i(n_i); + vector LogProb1_i(n_i); + // Linear predictor (pre-link) for positive component + vector P2_i(n_i); + // Response predictor (post-link) + // ObsModel = 0:3, 11:12: expected value of data, given that data is greater than zero -> E[D] = mu*phi + // ObsModel = 4 (ZANB): expected value ("mu") of neg-bin PRIOR to truncating Pr[D=0] -> E[D] = mu/(1-NB(0|mu,var))*phi ALSO Pr[D] = NB(D|mu,var)/(1-NB(0|mu,var))*phi + // ObsModel = 5 (ZINB): expected value of data for non-zero-inflation component -> E[D] = mu*phi + vector R2_i(n_i); + vector LogProb2_i(n_i); + vector maxJ_i(n_i); + vector diag_z(4); + matrix diag_iz(n_i,4); + diag_iz.setZero(); // Used to track diagnostics for Tweedie distribution (columns: 0=maxJ; 1=maxW; 2=lowerW; 3=upperW) + + // Likelihood contribution from observations + for(int i=0; i=0 & t_iz(i,z) 0 ){ + LogProb1_i(i) = log( R1_i(i) ); + }else{ + LogProb1_i(i) = log( 1-R1_i(i) ); + } + // Positive density likelihood -- models with continuous positive support + if( b_i(i) > 0 ){ // 1e-500 causes overflow on laptop + if(ObsModel(0)==0) LogProb2_i(i) = dnorm(b_i(i), R2_i(i), SigmaM(c_i(i),0), true); + if(ObsModel(0)==1) LogProb2_i(i) = dlnorm(b_i(i), log(R2_i(i))-pow(SigmaM(c_i(i),0),2)/2, SigmaM(c_i(i),0), true); // log-space + if(ObsModel(0)==2) LogProb2_i(i) = dgamma(b_i(i), 1/pow(SigmaM(c_i(i),0),2), R2_i(i)*pow(SigmaM(c_i(i),0),2), true); // shape = 1/CV^2, scale = mean*CV^2 + }else{ + LogProb2_i(i) = 0; + } + } + // Likelihood for Tweedie model with continuous positive support + if(ObsModel(0)==8){ + LogProb1_i(i) = 0; + //dPoisGam( Type x, Type shape, Type scale, Type intensity, Type &max_log_w_j, int maxsum=50, int minsum=1, int give_log=0 ) + LogProb2_i(i) = dPoisGam( b_i(i), SigmaM(c_i(i),0), R2_i(i), R1_i(i), diag_z, Options_vec(5), Options_vec(6), true ); + diag_iz.row(i) = diag_z; + } + if(ObsModel(0)==10){ + // Packaged code + LogProb1_i(i) = 0; + // dtweedie( Type y, Type mu, Type phi, Type p, int give_log=0 ) + // R1*R2 = mean + LogProb2_i(i) = dtweedie( b_i(i), R1_i(i)*R2_i(i), R1_i(i), invlogit(SigmaM(c_i(i),0))+1.0, true ); + } + // Likelihood for models with discrete support + if(ObsModel(0)==4 | ObsModel(0)==5 | ObsModel(0)==6 | ObsModel(0)==7 | ObsModel(0)==9){ + if(ObsModel(0)==5){ + // Zero-inflated negative binomial (not numerically stable!) + var_i(i) = R2_i(i)*(1.0+SigmaM(c_i(i),0)) + pow(R2_i(i),2.0)*SigmaM(c_i(i),1); + if( b_i(i)==0 ){ + LogProb2_i(i) = log( (1-R1_i(i)) + dnbinom2(Type(0.0), R2_i(i), var_i(i), false)*R1_i(i) ); // Pr[X=0] = 1-phi + NB(X=0)*phi + }else{ + LogProb2_i(i) = dnbinom2(b_i(i), R2_i(i), var_i(i), true) + log(R1_i(i)); // Pr[X=x] = NB(X=x)*phi + } + } + if(ObsModel(0)==6){ + // Conway-Maxwell-Poisson + LogProb2_i(i) = dCMP(b_i(i), R2_i(i), exp(P1_i(i)), true, Options_vec(5)); + } + if(ObsModel(0)==7){ + // Zero-inflated Poisson + if( b_i(i)==0 ){ + LogProb2_i(i) = log( (1-R1_i(i)) + dpois(Type(0.0), R2_i(i), false)*R1_i(i) ); // Pr[X=0] = 1-phi + Pois(X=0)*phi + }else{ + LogProb2_i(i) = dpois(b_i(i), R2_i(i), true) + log(R1_i(i)); // Pr[X=x] = Pois(X=x)*phi + } + } + if(ObsModel(0)==9){ + // Binned Poisson (for REEF data: 0=none; 1=1; 2=2-10; 3=>11) + /// Doesn't appear stable given spatial or spatio-temporal variation + vector logdBinPois(4); + logdBinPois(0) = logspace_add( log(1-R1_i(i)), dpois(Type(0.0), R2_i(i), true) + log(R1_i(i)) ); // Pr[X=0] = 1-phi + Pois(X=0)*phi + logdBinPois(1) = dpois(Type(1.0), R2_i(i), true) + log(R1_i(i)); // Pr[X | X>0] = Pois(X)*phi + logdBinPois(2) = dpois(Type(2.0), R2_i(i), true) + log(R1_i(i)); // SUM_J( Pr[X|X>0] ) = phi * SUM_J( Pois(J) ) + for(int j=3; j<=10; j++){ + logdBinPois(2) += logspace_add( logdBinPois(2), dpois(Type(j), R2_i(i), true) + log(R1_i(i)) ); + } + logdBinPois(3) = logspace_sub( log(Type(1.0)), logdBinPois(0) ); + logdBinPois(3) = logspace_sub( logdBinPois(3), logdBinPois(1) ); + logdBinPois(3) = logspace_sub( logdBinPois(3), logdBinPois(2) ); + if( b_i(i)==0 ) LogProb2_i(i) = logdBinPois(0); + if( b_i(i)==1 ) LogProb2_i(i) = logdBinPois(1); + if( b_i(i)==2 ) LogProb2_i(i) = logdBinPois(2); + if( b_i(i)==3 ) LogProb2_i(i) = logdBinPois(3); + } + LogProb1_i(i) = 0; + } + } + REPORT( diag_iz ); + + // Joint likelihood + jnll_comp(10) = -1 * (LogProb1_i * (Type(1.0)-PredTF_i)).sum(); + jnll_comp(11) = -1 * (LogProb2_i * (Type(1.0)-PredTF_i)).sum(); + jnll = jnll_comp.sum(); + Type pred_jnll = -1 * ( LogProb1_i*PredTF_i + LogProb2_i*PredTF_i ).sum(); + REPORT( pred_jnll ); + + //////////////////////// + // Calculate outputs + //////////////////////// + + // Number of output-years + int n_y = t_yz.col(0).size(); + + // Predictive distribution -- ObsModel(4) isn't implemented (it had a bug previously) + Type a_average = a_i.sum()/a_i.size(); + array P1_xcy(n_x, n_c, n_y); + array R1_xcy(n_x, n_c, n_y); + array P2_xcy(n_x, n_c, n_y); + array R2_xcy(n_x, n_c, n_y); + array D_xcy(n_x, n_c, n_y); + for(int c=0; c=0 & t_yz(y,z) Index_xcyl(n_x, n_c, n_y, n_l); + array Index_cyl(n_c, n_y, n_l); + array ln_Index_cyl(n_c, n_y, n_l); + Index_cyl.setZero(); + for(int c=0; c mean_Z_cym(n_c, n_y, n_m); + if( Options(2)==1 ){ + mean_Z_cym.setZero(); + int report_summary_TF = false; + for(int c=0; c mean_D_cyl(n_c, n_y, n_l); + array log_mean_D_cyl(n_c, n_y, n_l); + mean_D_cyl.setZero(); + for(int c=0; c effective_area_cyl(n_c, n_y, n_l); + array log_effective_area_cyl(n_c, n_y, n_l); + effective_area_cyl = Index_cyl / (mean_D_cyl/1000); // Correct for different units of Index and density + log_effective_area_cyl = log( effective_area_cyl ); + REPORT( effective_area_cyl ); + ADREPORT( effective_area_cyl ); + ADREPORT( log_effective_area_cyl ); + } + + // Reporting and standard-errors for covariance and correlation matrices + if( Options(5)==1 ){ + if( FieldConfig(0)>0 ){ + matrix L1_omega_cf = loadings_matrix( L_omega1_z, n_c, FieldConfig(0) ); + matrix lowercov_uppercor_omega1 = L1_omega_cf * L1_omega_cf.transpose(); + lowercov_uppercor_omega1 = convert_upper_cov_to_cor( lowercov_uppercor_omega1 ); + REPORT( lowercov_uppercor_omega1 ); + ADREPORT( lowercov_uppercor_omega1 ); + } + if( FieldConfig(1)>0 ){ + matrix L1_epsilon_cf = loadings_matrix( L_epsilon1_z, n_c, FieldConfig(1) ); + matrix lowercov_uppercor_epsilon1 = L1_epsilon_cf * L1_epsilon_cf.transpose(); + lowercov_uppercor_epsilon1 = convert_upper_cov_to_cor( lowercov_uppercor_epsilon1 ); + REPORT( lowercov_uppercor_epsilon1 ); + ADREPORT( lowercov_uppercor_epsilon1 ); + } + if( FieldConfig(2)>0 ){ + matrix L2_omega_cf = loadings_matrix( L_omega2_z, n_c, FieldConfig(2) ); + matrix lowercov_uppercor_omega2 = L2_omega_cf * L2_omega_cf.transpose(); + lowercov_uppercor_omega2 = convert_upper_cov_to_cor( lowercov_uppercor_omega2 ); + REPORT( lowercov_uppercor_omega2 ); + ADREPORT( lowercov_uppercor_omega2 ); + } + if( FieldConfig(3)>0 ){ + matrix L2_epsilon_cf = loadings_matrix( L_epsilon2_z, n_c, FieldConfig(3) ); + matrix lowercov_uppercor_epsilon2 = L2_epsilon_cf * L2_epsilon_cf.transpose(); + lowercov_uppercor_epsilon2 = convert_upper_cov_to_cor( lowercov_uppercor_epsilon2 ); + REPORT( lowercov_uppercor_epsilon2 ); + ADREPORT( lowercov_uppercor_epsilon2 ); + } + } + + // Synchrony + if( Options(6)==1 ){ + int n_z = yearbounds_zz.col(0).size(); + // Density ("D") or area-expanded total biomass ("B") for each category (use B when summing across sites) + matrix D_xy( n_x, n_y ); + matrix B_cy( n_c, n_y ); + vector B_y( n_y ); + D_xy.setZero(); + B_cy.setZero(); + B_y.setZero(); + // Sample variance in category-specific density ("D") and biomass ("B") + array varD_xcz( n_x, n_c, n_z ); + array varD_xz( n_x, n_z ); + array varB_cz( n_c, n_z ); + vector varB_z( n_z ); + vector varB_xbar_z( n_z ); + vector varB_cbar_z( n_z ); + vector ln_varB_z( n_z ); + vector ln_varB_xbar_z( n_z ); + vector ln_varB_cbar_z( n_z ); + array maxsdD_xz( n_x, n_z ); + array maxsdB_cz( n_c, n_z ); + vector maxsdB_z( n_z ); + varD_xcz.setZero(); + varD_xz.setZero(); + varB_cz.setZero(); + varB_z.setZero(); + varB_xbar_z.setZero(); + varB_cbar_z.setZero(); + maxsdD_xz.setZero(); + maxsdB_cz.setZero(); + maxsdB_z.setZero(); + // Proportion of total biomass ("P") for each location or each category + matrix propB_xz( n_x, n_z ); + matrix propB_cz( n_c, n_z ); + propB_xz.setZero(); + propB_cz.setZero(); + // Synchrony indices + matrix phi_xz( n_x, n_z ); + matrix phi_cz( n_c, n_z ); + vector phi_xbar_z( n_z ); + vector phi_cbar_z( n_z ); + vector phi_z( n_z ); + phi_xbar_z.setZero(); + phi_cbar_z.setZero(); + phi_z.setZero(); + // Calculate total biomass for different categories + for( int y=0; y CovHat( n_c, n_c ); + matrix CovHat( n_c, n_c ); + CovHat.setIdentity(); + CovHat *= pow(0.0001, 2); + if( FieldConfig(1)>0 ) CovHat += loadings_matrix(L_epsilon1_z, n_c, FieldConfig(1)) * loadings_matrix(L_epsilon1_z, n_c, FieldConfig(1)).transpose(); + if( FieldConfig(3)>0 ) CovHat += loadings_matrix(L_epsilon2_z, n_c, FieldConfig(3)) * loadings_matrix(L_epsilon2_z, n_c, FieldConfig(3)).transpose(); + // Coherence ranges from 0 (all factors are equal) to 1 (first factor explains all variance) + SelfAdjointEigenSolver > es(CovHat); + vector eigenvalues_c = es.eigenvalues(); // Ranked from lowest to highest for some reason + Type psi = 0; + for(int c=0; c diag_CovHat( n_c ); + vector log_diag_CovHat( n_c ); + for(int c=0; c