Catch Model Exploration

library(brms)

# load data
load("/Users/michaelgarbo/Documents/CATCH_SOURCES_CHAPTER/CATCH_all_working_R/data_objects/catch_data.rda")
# scaling and centring trip days
dh$zdays<-(dh$tripdays-mean(dh$tripdays))/sd(dh$tripdays)

# load model fit with brm
#hrXMS<- brm(data = dh, family = poisson,
#                c ~ 0 + intercept+xF+mF+sF+mF:sF+zdays+(1+zdays|vill)+(1|ht),
#                prior = c(prior(normal(0, 10), class = b, coef = intercept),
#                          prior(normal(0, 1), class = b),
#                          prior(cauchy(0, 1), class = sd)),
#                         iter=4000,cores=2,control=list(adapt_delta=0.9))
hrXMS<-readRDS(file="/Users/michaelgarbo/Documents/CATCH_SOURCES_CHAPTER/CATCH_all_working_R/brms_models/hr3XMiSz.rds")

I want to look at covariates that might give insight into ways catch data might become biased. The response variable (c) is the number of animals caught on a hunters’ most recent trip

Covariates are:

xF: whether hunter was accompanied by an extra person on the trip;0=no,1=yes

sF: season (early.dry, late.dry, rainy)

mF: hunting method used on trip (gun, trapper, both)

zdays: duration of trip in days scaled (subract mean, divide by SD)

The random effects are

vill: village ID

ht: hunter ID

The code to fit in package brms was this:

# not run
#hrXMS<- brm(data = dh, family = poisson,
#                c ~ 0 + intercept+xF+mF+sF+mF:sF+zdays+(1+zdays|vill)+(1|ht),
#                prior = c(prior(normal(0, 10), class = b, coef = intercept),
#                          prior(normal(0, 1), class = b),
#                          prior(cauchy(0, 1), class = sd)),
#                         iter=4000,cores=2,control=list(adapt_delta=0.9))
summary(hrXMS)

 Family: poisson 
  Links: mu = log 
Formula: c ~ 0 + intercept + xF + mF + sF + mF:sF + zdays + (1 + zdays | vill) + (1 | ht) 
   Data: dh (Number of observations: 252) 
Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup samples = 8000

Group-Level Effects: 
~ht (Number of levels: 208) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.42      0.04     0.34     0.50       3742 1.00

~vill (Number of levels: 18) 
                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)            0.53      0.13     0.32     0.82       4150 1.00
sd(zdays)                0.20      0.08     0.09     0.38       4464 1.00
cor(Intercept,zdays)    -0.10      0.42    -0.81     0.73       6436 1.00

Population-Level Effects: 
                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
intercept                2.07      0.17     1.74     2.40       2810 1.00
xF1                      0.29      0.08     0.14     0.44       9174 1.00
mFgun                   -0.22      0.14    -0.50     0.06       5749 1.00
mFtrapper               -0.16      0.15    -0.45     0.14       6038 1.00
sFlate.dry              -0.23      0.10    -0.43    -0.03       6921 1.00
sFrainy                  0.09      0.14    -0.18     0.36       8388 1.00
zdays                    0.28      0.08     0.12     0.43       6265 1.00
mFgun:sFlate.dry        -0.12      0.16    -0.44     0.20       5790 1.00
mFtrapper:sFlate.dry    -0.21      0.21    -0.62     0.20       7006 1.00
mFgun:sFrainy           -0.22      0.21    -0.61     0.19       6967 1.00
mFtrapper:sFrainy       -0.65      0.29    -1.23    -0.08       8185 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

Using leave one out cross validation and pareto-smoothed importance sampling (Vehtari, Gelman, and Gabry, 2017) shows that many data points influence the posterior estimates (i.e. acting like outliers).

loo(hrXMS)

Found 83 observations with a pareto_k > 0.7 in model 'hrXMS'. With this many problematic observations, it may be more appropriate to use 'kfold' with argument 'K = 10' to perform 10-fold cross-validation rather than LOO.

Computed from 8000 by 252 log-likelihood matrix

------
Monte Carlo SE of elpd_loo is NA.

Pareto k diagnostic values:
                         Count Pct.    Min. n_eff
(-Inf, 0.5]   (good)     73    29.0%   879       
 (0.5, 0.7]   (ok)       96    38.1%   206       
   (0.7, 1]   (bad)      75    29.8%   35        
   (1, Inf)   (very bad)  8     3.2%   7         
See help('pareto-k-diagnostic') for details.

I guess this is something to do with high variability between hunters and the fact that most hunters contributed only 1 data point. A model without hunter ID as a varying intercept has fewer problematic observations, but lower out of sample performance based on kfold CV.

Extra notes on loo_compare: elpd = estimated log of loo predictive density, unbiased estimate of log posterior predictive density for new data (i.e. if you have a new data point, log ppd of its value according to the model fit without it) p_loo = can be considered as effective number of parameters lpd - elpd; provides a measure of model complexity, how much the model was able to fit the data.

# hrXMS_b<-update(hrXMS,.~.-(1|ht))
loo(hrXMS_b)

Found 7 observations with a pareto_k > 0.7 in model 'hrXMS_b'. It is recommended to set 'reloo = TRUE' in order to calculate the ELPD without the assumption that these observations are negligible. This will refit the model 7 times to compute the ELPDs for the problematic observations directly.

Computed from 8000 by 252 log-likelihood matrix

------
Monte Carlo SE of elpd_loo is NA.

Pareto k diagnostic values:
                         Count Pct.    Min. n_eff
(-Inf, 0.5]   (good)     233   92.5%   370       
 (0.5, 0.7]   (ok)        12    4.8%   145       
   (0.7, 1]   (bad)        4    1.6%   54        
   (1, Inf)   (very bad)   3    1.2%   5         
See help('pareto-k-diagnostic') for details.

Comparing models with and without the covariates with kfold CV shows they all perform similarly in terms of out-of-sample prediction (the loo_compare output shows elpd_diff is less than 2* se_elpd_diff). I’m interpretting this as showing there is little value in knowing about covariates because hunters and villages themselves are so variable. I’m interpretting the posterior distributions as providing some evidence that catch is higher if hunters use both snares and guns, if they are accompanied by someone, and that trappers tend to have lower catch per day in the rainy season. Expected effects ignoring hunter and village-level variance look like this (conditional on base factor levels and mean trip duration, 3.9 days)) :

mz<-marginal_effects(hrXMS)
mzp<-plot(mz,ask=F,rug=T)

Including the effects of hunter and village variability reveals the large predictive uncertainty and better captures the spread of the data (points are raw data at the base factor levels) :

mz_re<-marginal_effects(hrXMS,re_formula=NULL)
mzp_re<-plot(mz_re,ask=F,rug=T,points=T,point_args=list(alpha=0.5,width=0.2))

An obvious problem with the model is that the relationship between catch and trip duration is crazy for longer trips, but there weren’t many trips longer than 10 days so I don’t know if I need to worry about this. In general the predictions for the range where most of the data fall look reasonable to me (see below).

Another (potentially more serious?) flaw is that it can predict a negative coeficient for trip duration - which obviously makes no sense. Here is the predicted relationships between catch and trip days for ‘new’ hunters for trips of 1 to 10 days, there are one or two lines that suggest a negative relationship, but again I don’t know if this is a problem that I really need to address given that I don’t care too much about making predictions or trying to get the most accurate parameter estimates, what I’m more interested in is basically that there is maybe some effects of covariates, and lots of variation across hunters and villages.

newd<-data.frame(sF='early.dry',mF=rep(c('both','gun','trapper'),2),xF=rep(c("0","1"),3))
pd<-make_conditions(newd,vars=c('mF','sF','xF'))
md<-marginal_effects(hrXMS, effects = 'zdays',conditions=pd,spaghetti=T,nsamples=200)
md_re<-marginal_effects(hrXMS, effects = 'zdays',conditions=pd,spaghetti=T,nsamples=200,re_formula=NULL)
ppd<-plot(md,points=T,plot=F)[[1]]
ppd_re<-plot(md_re,points=T,plot=F)[[1]]
pd1<-ppd+ggplot2::xlim(-1,2.04)+ggplot2::ylim(0,50)+ggplot2::ggtitle("ignore random effects")
pd2<-ppd_re+ggplot2::xlim(-1,2.04)+ggplot2::ylim(0,50)+ggplot2::ggtitle("include random effects")
pd1;pd2

residuals vs covariates checks look ok as far as I can tell, but there are probably more checks I should do

p0<-pp_check(hrXMS)

Using 10 posterior samples for ppc type 'dens_overlay' by default.

p02<-pp_check(hrXMS,x='zdays',type="error_scatter_avg_vs_x")

Using all posterior samples for ppc type 'error_scatter_avg_vs_x' by default.

p1a<-pp_check(hrXMS,type="scatter_avg_grouped",group="mF",nsamples=100)
p1b<-pp_check(hrXMS,type="scatter_avg_grouped",group="sF",nsamples=100)
p1c<-pp_check(hrXMS,type="scatter_avg_grouped",group="xF",nsamples=100)
p0;p02;p1a;p1b;p1c

Improving the model parameterisation: 1. Use prior predictive simulation to specify more appropriate priors 2. Compare models without random slopes for villages, see if predictions about zdays parameters is improved.

1. Better priors: using the argument sample_prior=‘only’ and looking at parameter estimates and implied predictions from the prior distributions to get an idea of priors that are more sensible and don’t have huge probability mass on estimates that are unrealistically high.

Using the same model structure with random slopes and intercepts but more regularising priors made for slightly less crazy prior estimates of hunters catch, but didn’t solve the problem that the relationship with trip duration could be predicted as negative for some hunters.

plot(pmp,ask=F,points=T,point_args=list(width=0.3,alpha=0.6))[[4]]+ggplot2::xlim(-1,2)+ggplot2::ylim(0,40000) # zooming in on the zdays the points are for baseline factor levels only

# including random effects priors can allow negative coefficients for tripdays
rmps<-marginal_effects(m2,re_formula=NULL,spaghetti=T,nsamples=500,plot=F)
plot(rmps,ask=F,points=T,point_args=list(width=0.3,alpha=0.6))[[4]]+ggplot2::xlim(-1,2.04)+ggplot2::ylim(0,400)

# fit the model with better priors and look at implied predictions of the posterior
hrXMS3<- brm(data = dh, family = poisson,
                  c ~ 0 + intercept +xF+mF+sF+mF:sF+zdays+(1+zdays|vill)+(1|ht),
                  prior = c(prior(normal(0, 5), class = b, coef = intercept),
                            prior(normal(0, 0.5), class = b,coef='zdays'),
                            prior(normal(0, 0.25), class = b),
                            prior(exponential(2), class = sd)),iter=4000,cores=2,control=list(adapt_delta=0.9),file="/Users/michaelgarbo/Documents/CATCH_SOURCES_CHAPTER/CATCH_all_working_R/brms_models/hrXMiS3")

m2m<-hrXMS3

# spaghetti plot shows one or two zdays coefficients are still negative and there are still some crazy high predictions for longer trips
rmpsP<-marginal_effects(m2m,re_formula=NULL,spaghetti=T,nsamples=500,plot=F)
plot(rmpsP,ask=F,points=T,point_args=list(width=0.3,alpha=0.6))[[4]]+ggplot2::xlim(-1,2.04)+ggplot2::ylim(0,200)

2. Changing the model structure to leave out a varying slope at the village level (i.e. the relationship between catch and trip duration is now assumed to be the same across villages) produced only positive predicted relationships and controlled the high values of catch estimated for longer trips. Worth noting that the posterior distribution of the correlation parameter between intercept and slope for trip duration had had a high density over 0 and large variance - so basically the model was really not sure about what, if any, the correlation should be between the average catch at a given village, and the increase you’d expect with longer trips. Since there was also an intercept for each hunter (and many hunters only contribute one data point) it makes sense that there really isn’t enough data to include random slope for villages.

Some basic plots of the error distributions doesn’t flag up anything to me, although the predictive error is pretty large in general.

Next look at the sensitivity of predictions to the priors used, I fit another model with the same structure but less informative priors and visualised the prior distributions and posterior distributions of the parameter estimates. Both the models ended up with similar parameter estimates but the interaction term parameters seemed like they might have been a bit constrained by the priors, so I refit with a wider prior on the coefficients for the fixed effects, (and the predictions etc still looked sensible).

plots of posterior parameter distributions, pretty good evidence that trips that had an extra person had higher catch than those which didn’t (makes sense, especially if hunters report their helpers’ catch as part of their own)

hx

[[1]]

hx2

[[1]]

Also that hunting trips where both guns and trapping were used had higher catch than those only using one method, not shown here the estimate for trappers suggests no difference between trappers and gun hunters

hgun

[[1]]

hgun2

[[1]]

no evidence that trappers did more poorly in the rainy season than the early dry season (which is the factor baseline)

hrtraprainy

[[1]]

hrtraprainy2

[[1]]

and hunters catch more on longer trips (duh).

hdays

[[1]]

hdays2

[[1]]

Changing the priors to be less restrictive on the fixed effects parameters seems like a good idea so I try that

hrXMS3b<- brm(data = dh, family = poisson,
              c ~ 0 + intercept +xF+mF+sF+mF:sF+zdays+(1|vill)+(1|ht),
              prior = c(prior(normal(0, 5), class = b, coef = intercept),
                        prior(normal(0, 0.5), class = b),
                        prior(exponential(2), class = sd)),
              iter=4000,cores=2,control=list(adapt_delta=0.9),sample_prior=T,file="/Users/michaelgarbo/Documents/CATCH_SOURCES_CHAPTER/CATCH_all_working_R/brms_models/hrXMiS3b")

summary(hrXMS3b)

and I’m still quite surprised by the posterior of the interaction parameter for trappers in the rainy season, although the other posterior parameters have more mass on 0 but have large uncertainty. Then I learned that the estimate of an interaction has twice the standard error of an estimate for a main effect, so actually this is inevitable. Also it means that I almost certainly don’t have anywhere near enough data to estimate an interaction effect, although I guess that doesn’t mean I should not have it in the model.

looking back at the raw data shows there are only 8 data points for trappers in the rainy season, and their catch is quite variable and there is not obviously a big interaction effect but conditioning on the interaction might help get more reasonable parameter distributions for the main effects.

datap1<-ggplot(dhh,aes(x=mF,y=cpd,alpha=0.5,colour=sF))+geom_point(position=position_jitterdodge())+geom_boxplot(varwidth=T)
datap2<-ggplot(dhh,aes(x=sF,y=cpd,alpha=0.5,colour=mF))+geom_point(position=position_jitterdodge())+geom_boxplot(varwidth=T)
datap3<-ggplot(dhh,aes(x=mF,y=c,alpha=0.5,colour=sF))+geom_point(position=position_jitterdodge())+geom_boxplot(varwidth=T)
datap4<-ggplot(dhh,aes(x=sF,y=c,alpha=0.5,colour=mF))+geom_point(position=position_jitterdodge())+geom_boxplot(varwidth=T)
datap1
datap2

datap3

datap4

Given that there isn’t much data for this condition in the interaction term, I guess theres a risk that whatever priors I choose may influence the posterior somewhat. So I’m going to stick with the normal(0,0.5) option for now, and be cautious about interpreting the effects of the interaction between methods and season. Also I want to look at the out of sample performance for these models….

The LOO-PSIS check still shows that individual data points have substantial impacts on the posterior distributions, so I use kfold cross-validation to see whether taking out any of the covariates makes a big difference in terms of how model fit.

loo_compare(hrXMS3pc,hrXMS3b,hrXMS3a,criterion='kfold')

         elpd_diff se_diff
hrXMS3pc   0.0       0.0  
hrXMS3a   -2.7       9.1  
hrXMS3b  -16.7       9.6  

loo_compare(hrXMS3a,hrXMS4,hrMS5,hrMX4,criterion='kfold')

        elpd_diff se_diff
hrXMS3a   0.0       0.0  
hrXMS4   -2.7      11.7  
hrMX4   -12.6       8.7  
hrMS5   -19.2      25.2  

The models with different priors, and with and without the interaction term and dropping the fixed effects don’t seem to vary much in their out of sample predictive performance. So just sticking with the original model for now….

Start to explore the posterior in more detail, e.g. look at village and hunter level variation in intercepts

The overall marginal effects of this model (predicted catch for the baseline factor levels) look pretty similar to those in the first model which had a random slope and weaker priors.

mp<-marginal_effects(mod,re_formula=NA,plot=F)
pp<-plot(mp,ask=F,points=T,point_args=list(width=0.3,alpha=0.6))

pp[[4]]+ggplot2::xlim(-1,2)+ggplot2::ylim(0,40) # zooming in on the zdays note the points are for baseline factor levels only

# including the random effects and spaghetti plot 
mps<-marginal_effects(mod,re_formula=NULL,spaghetti=T,nsamples=500,plot=F)
plot(mps,ask=F,points=T,point_args=list(width=0.3,alpha=0.6))[[4]]+ggplot2::xlim(-1,2.04)+ggplot2::ylim(0,40)