Reputation: 67
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)
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)
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