# Multiplying glmer binomial predicted probabilities to account for overall proportions in population

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))