bribina
bribina

Reputation: 67

How to display confidence intervals of bootstrapping GLM results for each factor?

tl;dr: I bootstrapped a generalized linear model. I would like to see the confidence interval for each factor instead of the CIs of reference levels.

I ran a generalized linear model. My dependent variable (body-phosphorus content) is numeric and continuous, my factors are fixed. One factor is binary (diet) and the other is 3 categorical units (reproductive output) - 2 of the units are the same samples but different values. To be more specific about the categorical factor, I have snails that were not reproductive, snails that were reproductive without their brood mass, and snails that were reproductive with their brood mass. The snails that were reproductive are the same snails but their values are different because of the inclusion or exclusion of the brood mass. I do not believe all results that I listed here are necessary to answer my question so I suggest that readers skip to the bootstrapping portion below. Code:

library("glmmTMB")
mod <- glmmTMB(body-phosphorus_content ~ diet*Reproductive_output, data = data, family = "Gamma")
summary(mod)
car::Anova((mod), type = c("III"))

Output:

> summary(mod)

Call:
glm(formula = Weight_Percentage ~ P_Treatment * Reproductive_output, 
    family = "Gamma", data = kdata)

Coefficients:
                                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)                             1.35516    0.08267  16.393  < 2e-16 ***
P_TreatmentLow                          0.23227    0.11702   1.985   0.0486 *  
Reproductive_outputBwoE                -0.03811    0.10862  -0.351   0.7261    
Reproductive_outputBwE                 -0.78160    0.08832  -8.850 5.31e-16 ***
P_TreatmentLow:Reproductive_outputBwoE -0.11308    0.16396  -0.690   0.4912    
P_TreatmentLow:Reproductive_outputBwE  -0.22673    0.12712  -1.784   0.0760 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.1116333)

    Null deviance: 64.223  on 200  degrees of freedom
Residual deviance: 25.058  on 195  degrees of freedom
AIC: 116.7

Number of Fisher Scoring iterations: 5

> car::Anova((mod), type = c("III"))
Analysis of Deviance Table (Type III tests)

Response: Weight_Percentage
                                LR Chisq Df Pr(>Chisq)    
P_Treatment                        3.911  1    0.04796 *  
Reproductive_output              164.942  2    < 2e-16 ***
P_Treatment:Reproductive_output    3.571  2    0.16770    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

My model showed that my one factor had unequal variances across the variants of Reproductive_output. My code:

library("DHARMa")
sim <- simulateResiduals(mod)
plot(sim)
plotResiduals(sim, data$Reproductive_output) plotResiduals(sim, data$diet)

Residual plot for Gamma Gamma Levene test for categorical factor Gamma Levene test for binary factor

I bootstrapped my model - my code:

library("car")
betahat.bootG <- Boot(mod, R=1000) 
summary(betahat.bootG) # default summary
confint(betahat.bootG) 
hist(betahat.bootG)

My output:

> summary(betahat.bootG) # default summary

Number of bootstrap replications R = 1000 
                                        original    bootBias   bootSE  bootMed
(Intercept)                             1.355159  0.00464240 0.096756  1.35591
P_TreatmentLow                          0.232265  0.00042426 0.125990  0.23337
Reproductive_outputBwoE                -0.038105 -0.00262016 0.114032 -0.04163
Reproductive_outputBwE                 -0.781598 -0.00126816 0.105097 -0.77561
P_TreatmentLow:Reproductive_outputBwoE -0.113083 -0.00139765 0.156818 -0.11464
P_TreatmentLow:Reproductive_outputBwE  -0.226726 -0.00423723 0.137830 -0.23268
> confint(betahat.bootG) 

Bootstrap bca confidence intervals
                                               2.5 %     97.5 %
(Intercept)                             1.1677409198  1.5494715
P_TreatmentLow                         -0.0008315274  0.4716950
Reproductive_outputBwoE                -0.2462003152  0.1973684
Reproductive_outputBwE                 -1.0149864479 -0.5963646
P_TreatmentLow:Reproductive_outputBwoE -0.4154117069  0.1894941
P_TreatmentLow:Reproductive_outputBwE  -0.4856801170  0.0716654

> hist(betahat.bootG)

Histogram of bootstrapping models

To make sure, I understand the output, I will not be able to see the confidence intervals for each of the variations (for example: Reproductive_outputBwoE, Reproductive_outputBwE, and Reproductive_outputNon-brooding), right? What code could I use to potentially see just the factors similar to the layout of my ANOVA table?

Thank you!

Upvotes: 0

Views: 50

Answers (0)

Related Questions