2
$\begingroup$

I am looking at the proportion of fishing trips made using a particular gear (LE_MSZ2) in two management phases (FT_PHASE) to look for change. Each observation is a number of fishing trips made using a gear in a phase. I have fitted a binomial glmer to the data, however the predicted proportions of gear used are the proportions IF a particular gear is used. So the predicted proportions in each management phase are not summing to 1 as expected. This allows me to answer the question: if a boat uses gear A does the proportion of trips it uses gear A in increase or decrease with the change in management phase. However, it doesn't let me answer overall for all boats has the proportion of trips using gear A increased or decreased with the change in management phase. Looking at the plots it is slightly confusing because the increases observed for some gears appear to be much larger than the decreases observed elsewhere in other gears. This is because the model does not appear to be taking into account the proportion of trips using each gear in each phase by all vessels.

In the plots below I have simply multiplied the model predictions by the proportion of trips made in each gear and each management phase to generate the results I was expecting. But is there a way to do this within the model - via weighting perhaps? OR can I have a two step model that first calculates the proportion of trips using each gear in each phase and then calculates the proportion of total trips using each gear when a gear is chosen (the second is the model I currently believe I have).

        summary(model)
        Generalized linear mixed model fit by maximum likelihood (Laplace
          Approximation) [glmerMod]
         Family: binomial  ( logit )
        Formula: cbind(FT_REF, other) ~ (1 | VE_NAME) + (1 | ID) + FT_PHASE *  
            LE_MSZ2
           Data: counts_APr
        Control: 
        glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 200000))
        
             AIC      BIC   logLik deviance df.resid 
           886.1    919.2   -431.1    862.1      104 
        
        Scaled residuals: 
             Min       1Q   Median       3Q      Max 
        -0.50265 -0.07522  0.01066  0.11327  0.80610 
        
        Random effects:
         Groups  Name        Variance           Std.Dev.    
         ID      (Intercept) 8.6553220865733902 2.9419928767
         VE_NAME (Intercept) 0.0000000000008947 0.0000009459
        Number of obs: 116, groups:  ID, 116; VE_NAME, 30
    
    
    Fixed effects:
                                       Estimate Std. Error z value   Pr(>|z|)    
    (Intercept)                         -3.3326     1.2419  -2.683    0.00729 ** 
    FT_PHASEAuditable                    0.3284     3.2225   0.102    0.91883    
    LE_MSZ2(90,100]                     -0.2653     2.1478  -0.124    0.90171    
    LE_MSZ2(100,120]                     6.1069     1.3753   4.441 0.00000897 ***
    LE_MSZ2(120,125]                     2.9676     1.5308   1.939    0.05255 .  
    LE_MSZ2(125,140]                     0.7134     1.4928   0.478    0.63270    
    FT_PHASEAuditable:LE_MSZ2(90,100]   -1.6123     4.3029  -0.375    0.70788    
    FT_PHASEAuditable:LE_MSZ2(100,120]  -0.4301     3.3308  -0.129    0.89726    
    FT_PHASEAuditable:LE_MSZ2(120,125]   0.9572     3.4889   0.274    0.78382    
    FT_PHASEAuditable:LE_MSZ2(125,140]   0.7027     3.4339   0.205    0.83786    
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library(ggeffects)

mydf <- ggpredict(model, terms = c("FT_PHASE", "LE_MSZ2"))

ggplot(data = mydf, aes(x=group, y=predicted, color = x)) +
  geom_point(size = 3, position=position_dodge(width=0.5)) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    width = 0.3,
    linetype = "solid",
    position=position_dodge(width=0.5)) +
  labs(x="\nMesh size (mm)",y="\nProportion of trips",color="Phase",title="glmer model i.e. predicted proportion \n when vessel uses mesh size category") +
  theme_bw()  +
  theme(axis.text = element_text(size = 10)) 

counts_APr <- counts_APr %>% 
  mutate(
    PHASE_SUM = case_when(
      FT_PHASE == "Pre-scheme" ~ "3585",
      FT_PHASE != "Pre-scheme" ~ "2504"
    )
  )
# calculate the overall proportions of gear use in each phase
multiplier <- counts_APr %>% 
  group_by(FT_PHASE,LE_MSZ2) %>% 
  summarise(trips = sum(FT_REF),
            sum = as.numeric(PHASE_SUM[1]),
            prop = trips/sum)
left  <- merge(mydf, multiplier, by.x = c("x","group"),by.y = c("FT_PHASE","LE_MSZ2"), all.x = TRUE)

# merge proportions with model predictions
mydf <- left %>% 
  mutate(predicted = predicted*prop,
         conf.low = conf.low*prop,
         conf.high = conf.high*prop)

# plot expected data
ggplot(data = mydf, aes(x=group, y=predicted, color = x)) +
  geom_point(size = 3, position=position_dodge(width=0.5)) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    width = 0.3,
    linetype = "solid",
    position=position_dodge(width=0.5)) +
  labs(x="\nMesh size (mm)",y="\nProportion of trips",color="Phase",title = "adjusted for proportion of sample trips") +
  theme_bw()  +
  theme(axis.text = element_text(size = 10)) 

This is the model predicted output - where you can see the increases in 120,125 and 125, 140 not seem to be matched by a decrease anywhere

This is the adjusted model predicted output which looks more as expected

$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.