Figs Revisited
Predictor FIG.DBH_m
Response Bird \(\rightarrow\) count data \(\rightarrow\)
Poisson family
look
figs <- read.csv("fig_visitors.csv")
#unit confirmation
summary(figs$FIG.DBH)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 54.2 82.3 119.1 143.8 140.1 568.8range(figs$FIG.DBH, na.rm = T)## [1] 54.2 568.8#54.2 to 568.8 suggest the data are in cm...not m. So let's convert for last question. The reason I noticed this is because the graphs of posterior draws had 2 humps.
figs$FIG.DBH_m <- figs$FIG.DBH/100
ggplot(figs, aes(FIG.DBH_m, Bird, color="pink"))+
geom_point()+
theme(legend.position = "none")+
labs(x="Fig tree diameter (m)", y="Bird visitors")
model
m <- brm(Bird ~ FIG.DBH_m, data=figs, family="poisson", refresh=0)## Running "C:/PROGRA~1/R/R-45~1.2/bin/x64/Rcmd.exe" SHLIB foo.c
## using C compiler: 'gcc.exe (GCC) 14.3.0'
## gcc -I"C:/PROGRA~1/R/R-45~1.2/include" -DNDEBUG -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/Rcpp/include/" -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/" -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/unsupported" -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/BH/include" -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/src/" -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/" -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppParallel/include/" -DRCPP_PARALLEL_USE_TBB=1 -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include "C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp" -std=c++1y -I"C:/rtools45/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -std=gnu2x -mfpmath=sse -msse2 -mstackrealign -c foo.c -o foo.o
## cc1.exe: warning: command-line option '-std=c++14' is valid for C++/ObjC++ but not for C
## In file included from C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Core:19,
## from C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Dense:1,
## from C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [C:/PROGRA~1/R/R-45~1.2/etc/x64/Makeconf:289: foo.o] Error 1posterior draws
posterior_draws <- as.data.frame(m)
head(posterior_draws)## b_Intercept b_FIG.DBH_m Intercept lprior lp__
## 1 2.201474 0.2719160 2.592415 -1.918090 -121.3381
## 2 2.298113 0.2522930 2.660841 -1.919937 -121.1319
## 3 2.298113 0.2522930 2.660841 -1.919937 -121.1319
## 4 2.307213 0.2004029 2.595338 -1.918149 -121.7838
## 5 2.292885 0.2007634 2.581528 -1.917888 -122.0093
## 6 2.272656 0.2757753 2.669145 -1.920229 -121.8948plot posterior draws
plot(density(posterior_draws$b_Intercept), main="Intercept's posterior draws", col="tomato")
plot(density(posterior_draws$b_FIG.DBH_m), "Slope's posterior draws", col="steelblue")
summarize posteriors
mean_intercept <- mean(posterior_draws$b_Intercept)
quantile(posterior_draws$b_Intercept, c(0.025, 0.5, 0.975))## 2.5% 50% 97.5%
## 2.110271 2.283895 2.449672mean_slope <- mean(posterior_draws$b_FIG.DBH_m)
quantile(posterior_draws$b_FIG.DBH_m, c(0.025, 0.5, 0.975))## 2.5% 50% 97.5%
## 0.1672194 0.2386419 0.3083459scale slope & summarize
Convert slope posterior to interpretative scale:
scaled_slope <- exp(posterior_draws$b_FIG.DBH_m)
plot(density(scaled_slope), main="Multiplicative effect per +1 DBH", xlab="exp(posterior slopes (FIG.DBH_m))")
print(c("Mean scaled slope:", mean(scaled_slope)))## [1] "Mean scaled slope:" "1.26924384761423"quantile(scaled_slope, c(0.025, 0.5, 0.975))## 2.5% 50% 97.5%
## 1.182014 1.269524 1.361172The quantile(scaled_slope) values are not percent change
yet. Those are computations for multiplicative effects. Now we have to
translate them.
Technically, each of the above values means:
For every +1 m of increased
FIG.DBH_mthere is a multiplicative increase between: 1.2 - 1.4 (95% credible interval)Birdcount, median = 2.3.
About 1.27 times more birds per unit fig tree diameter.
Convert to %
percents <- (scaled_slope-1)*100
quantile(percents, c(0.025, 0.5, 0.975))## 2.5% 50% 97.5%
## 18.20136 26.95239 36.11717For every 1 m increase in DBH, the expected number of bird visitors increases by a median of 26.8% with a 95% credible interval of 18.0 to 35.6%.
Note that Poisson regression gives you multiplicative slopes, not additive slopes like linear regression.
interval question
What is the prediction of bird visitors between 0.5 - 5 m DBH?
To answer this we compute the effect size by taking posterior draws of intercepts and slopes and subtracting two exponentiated predictions:
fixef(m)## Estimate Est.Error Q2.5 Q97.5
## Intercept 2.2845512 0.08610122 2.1102707 2.4496719
## FIG.DBH_m 0.2377688 0.03614015 0.1672194 0.3083459lower <- 0.5
upper <- 5
a <- 2.28 #from fixef
b <- 0.2375
lambda_0.5 <- exp(posterior_draws$b_Intercept+posterior_draws$b_FIG.DBH_m*lower)
lambda_5 <- exp(posterior_draws$b_Intercept+posterior_draws$b_FIG.DBH_m*upper)
effect_size <- lambda_5-lambda_0.5
quantile(effect_size, c(0.025, 0.5, 0.975))## 2.5% 50% 97.5%
## 13.36269 21.24769 30.56198Increasing fig DBH from 0.5 m to 5 m is associated with a median increase of 21 bird visits, with a 95% credible interval of 13 - 30 birds.
plot(density(effect_size), main="Expected increase in birds: DBH 0.5 - 5 m", xlab="Increase in expected bird visits")