Data analysis example: Hurricane lizards
Peter Ralph
8 February 2021 – Advanced Biological Statistics
Hurricane lizards
Data from: Donihue, C.M., Herrel, A., Fabre, AC. et al. Hurricane-induced selection on the morphology of an island lizard. Nature 560, 88–91 (2018).
Anolis scriptus morphology and performance from before and after hurricanes
Lizard morphology data collected from two visits to Pine Cay and Water Cay in
the Turks and Caicos Islands. Anolis scriptus lizards were surveyed on both
islands four days before hurricane Irma and again, six weeks later, after the
islands had been hit by Hurricanes Irma and Maria. We measured morphology and
performance and found significant differences in the "before" and "after"
hurricane lizard population morphology. All linear measurements are in MM,
surface area measurements are in MM2 and force in N. Counts correspond to the
number of lamellar scales on the forelimb and hind limb toepads.
Here’s the data: ../Datasets/Hurricane_lizards/. First we’ll read it in, reorder the levels of the Hurricane factor (so that “before” comes before “after”), change units on SVL to centimeters so it’s on the same scale as the other variables (useful below), and drop unused variables.
lizards <- read.csv("../Datasets/Hurricane_lizards/lizards.csv", header=TRUE, stringsAsFactors=TRUE)
lizards$Hurricane <- factor(lizards$Hurricane, levels=c("Before", "After"))
lizards$SVL_cm <- lizards$SVL / 10
lizards <- lizards[,setdiff(names(lizards), c(paste0("FingerArea", 1:3), paste0("ToeArea", 1:3), "SumFingers", "SumToes", "MaxFingerForce", "SVL"))]
head(lizards)## ID Hurricane Origin Sex Femur Tibia Metatarsal LongestToe Humerus Radius Metacarpal LongestFinger FingerCount ToeCount MeanFingerArea MeanToeArea SVL_cm
## 1 537 After Pine Cay Male 10.39 11.87 7.52 7.43 8.66 7.99 2.22 3.19 10 12 1.3327 2.433 4.869
## 2 539 After Pine Cay Female 8.66 9.79 6.18 6.20 8.01 6.51 2.38 3.55 10 13 0.9613 1.518 4.031
## 3 540 After Pine Cay Male 12.87 14.76 9.45 9.58 11.72 9.54 3.54 5.09 14 15 2.6313 4.098 5.830
## 4 541 After Pine Cay Female 8.55 10.29 6.60 6.26 7.43 6.60 2.79 3.55 11 12 1.1777 1.879 4.315
## 5 542 After Pine Cay Female 10.26 11.02 6.89 7.02 7.71 7.25 2.52 3.37 11 13 1.3843 2.530 4.551
## 6 543 After Pine Cay Female 10.02 10.78 6.85 7.18 8.45 7.15 2.39 3.26 12 14 1.4260 2.036 4.697Look at the data
First, let’s have a look at the sample sizes:
## , , Hurricane = Before
##
## Origin
## Sex Pine Cay Water Cay
## Female 15 18
## Male 18 20
##
## , , Hurricane = After
##
## Origin
## Sex Pine Cay Water Cay
## Female 17 22
## Male 29 25Our preliminary brainstorming suggested that we simply look at how the means (and distributions) of the variables changed. (For instance, did all the lizards longer than 50cm get blown away?) A good way to visualize this is with some boxplots, separated by before/after hurricane and island. To do this, we’ll first reshape the data frame, making a “long” version.
(ggplot(lizards_long) + geom_boxplot(aes(x=name, y=value, col=Hurricane))
+ theme(axis.text.x=element_text(angle=90))
+ facet_wrap(~ Origin)
)## Warning: Removed 4 rows containing non-finite values (stat_boxplot).
t-tests
There’s nothing super obvious going on there, but just to have a look at the differences, we could run some \(t\) tests (basically, one variable, comparing “before” to “after”). This ignores lots of information we have (e.g., island, sex, other variables…), but is a good first step.
tt_list <- lapply(varnames, function (varn) {
t.test(lizards[[varn]][lizards$Hurricane == "Before"],
lizards[[varn]][lizards$Hurricane == "After"])
} )
names(tt_list) <- varnames
t_results <- data.frame(variable = names(tt_list))
t_results$statistic <- t_results$p.value <- NA
for (j in 1:nrow(t_results)) {
stopifnot(names(tt_list)[j] == t_results$variable[j])
t_results$statistic[j] <- tt_list[[j]][["statistic"]]
t_results$p.value[j] <- tt_list[[j]][["p.value"]]
}
# adjust p-values for multiple testing
t_results$adj_p <- p.adjust(t_results$p.value)Here’s the results:
## variable p.value statistic adj_p
## 1 SVL_cm 0.50488 0.6686 1.0000
## 2 Femur 0.01555 2.4492 0.2022
## 3 Tibia 0.40788 0.8302 1.0000
## 4 Metatarsal 0.51235 0.6569 1.0000
## 5 LongestToe 0.04570 2.0177 0.5027
## 6 Humerus 0.47341 -0.7188 1.0000
## 7 Radius 0.59784 0.5287 1.0000
## 8 Metacarpal 0.69291 -0.3957 1.0000
## 9 LongestFinger 0.76508 0.2994 1.0000
## 10 FingerCount 0.02528 -2.2604 0.3033
## 11 ToeCount 0.52614 -0.6354 1.0000
## 12 MeanFingerArea 0.11483 -1.5863 1.0000
## 13 MeanToeArea 0.29694 -1.0468 1.0000Results of preliminary one-variable analysis
The overall means of these traits has not shifted a large amount as a result of the hurricane: examination of the boxplots (above) shows that the distributions, before and after, by island, largely or entirely overlap, and the mean values are at most shifted by a small amount (at most 6%, for femur length). None of these differences are statistically significant after adjusting for multiple corrections; however, this may be due to a lack of power (i.e., there may be a difference in trait values of around 5%, but the sample size is not sufficient to tell for sure).
Pairwise plots
This maybe should have been our first plot - pair plots often let us see a lot more structure, and this is no exception. Everything is correlated… with size! In class we put these plots in a pdf, to better zoom in. Below “black” is the first level in the title, and “red” is the second.
Think about the questions
Goals
- describe how traits differ before/after hurricanes
- account for confounding factor: overall size
- also possibly account for differences between islands and sexes
For instance, we’d like to say something like:
Remaining male lizards of average size on Pine Cay had toepads that were X% larger by area (95% CI of L-U%). This represents roughly Y standard deviations of toepad area. Toepad area did not differ substantially by island, or by sex after accounting for length.
Fit some models
lm
We’ll start off with lm( ) because (a) why not, and (b) to give us a consistency check for later.
##
## Call:
## lm(formula = MeanToeArea ~ Hurricane + SVL_cm, data = lizards)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6388 -0.1517 -0.0032 0.1464 0.8583
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.3819 0.1545 -28.37 < 2e-16 ***
## HurricaneAfter 0.2386 0.0384 6.22 4.2e-09 ***
## SVL_cm 1.3654 0.0304 44.93 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.243 on 160 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.927, Adjusted R-squared: 0.926
## F-statistic: 1.02e+03 on 2 and 160 DF, p-value: <2e-16plot(MeanToeArea ~ SVL_cm, data=lizards, pch=20, col=Hurricane)
abline(-4.14327, 1.3654, col='red')
abline(-4.14327-0.23863, 1.3654)
On average, lizards of the same length after the hurricane had toe pad areas 0.24 mm2 larger than before the hurricane (SE .04 mm2).
To build out the model better, and to more flexibly query it, we’ll switch to brms. We start simple for (a) a consistency check, and (b) if bad things happen with more complex models, we can step back and figure out what happened.
brm_lm <- brm(MeanToeArea ~ Hurricane + SVL_cm, data=lizards,
prior=set_prior("normal(0, 5)", class="b"))## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Start sampling## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: MeanToeArea ~ Hurricane + SVL_cm
## Data: lizards (Number of observations: 163)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -4.38 0.16 -4.70 -4.07 1.00 4548 2873
## HurricaneAfter 0.24 0.04 0.16 0.31 1.00 3882 2948
## SVL_cm 1.37 0.03 1.30 1.43 1.00 4647 3155
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.24 0.01 0.22 0.27 1.00 4229 3060
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Additive effects of other variables
We should consider possible effects of the other factors: what if larger toes were just due to higher mortality on the big-toed island?
brm_lm2 <- brm(MeanToeArea ~ Hurricane + SVL_cm + Sex + Origin,
data=lizards,
prior=set_prior("normal(0, 5)", class="b"))## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Start sampling## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: MeanToeArea ~ Hurricane + SVL_cm + Sex + Origin
## Data: lizards (Number of observations: 163)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -4.25 0.29 -4.82 -3.68 1.00 2696 2801
## HurricaneAfter 0.23 0.04 0.16 0.31 1.00 3991 3014
## SVL_cm 1.34 0.06 1.21 1.47 1.00 2685 2614
## SexMale 0.03 0.08 -0.12 0.19 1.00 2577 2397
## OriginWaterCay -0.03 0.04 -0.11 0.04 1.00 3653 2466
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.25 0.01 0.22 0.27 1.00 3935 2761
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(brm_lm2, effects="SVL_cm:Hurricane", conditions=data.frame(Origin=levels(lizards$Origin), cond__=levels(lizards$Origin)))
Varying slopes
If big toepads were important for big lizards but not so much for small ones, that’d change the slope, too. We might as well allow the slope to vary with the other variables, also? (Perhaps we’ll regret that…)
brm_lm3 <- brm(MeanToeArea ~ SVL_cm * (Hurricane + Sex + Origin),
data=lizards,
prior=set_prior("normal(0, 5)", class="b"))## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Start sampling## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: MeanToeArea ~ SVL_cm * (Hurricane + Sex + Origin)
## Data: lizards (Number of observations: 163)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -2.40 0.68 -3.74 -1.04 1.00 1495 2095
## SVL_cm 0.92 0.15 0.62 1.23 1.00 1517 2095
## HurricaneAfter -0.50 0.30 -1.08 0.07 1.00 2493 2345
## SexMale -1.99 0.74 -3.49 -0.53 1.00 1520 1896
## OriginWaterCay -0.06 0.31 -0.65 0.54 1.00 1831 2329
## SVL_cm:HurricaneAfter 0.15 0.06 0.04 0.27 1.00 2489 2481
## SVL_cm:SexMale 0.44 0.16 0.12 0.77 1.00 1505 1871
## SVL_cm:OriginWaterCay 0.00 0.06 -0.12 0.12 1.00 1843 2342
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.24 0.01 0.21 0.27 1.00 3470 2531
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Huh, where’s our hurricane effect gone?
conditional_effects(brm_lm3, effects="SVL_cm:Hurricane", conditions=data.frame(Origin=levels(lizards$Origin), cond__=levels(lizards$Origin)))
Don’t forget to think about what the reference categories are!!
condf <- expand.grid(Origin=levels(lizards$Origin), Sex=levels(lizards$Sex))
condf$cond__ <- paste(condf$Origin, condf$Sex)
plot(conditional_effects(brm_lm3, effects="SVL_cm:Hurricane", conditions=condf, plot=FALSE), facet_args=list(ncol=2), points=TRUE)
Going back to the data:
We have no power to distinguish an effect of sex from a nonlinear effect of SVL. 
Allowing a nonlinear term on SVL should be equivalent to including sex, and be less problematic. (Recall that I( ) in a formula means “no, really, this thing exactly”.)
brm_lm4 <- brm(MeanToeArea ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2),
data=lizards,
prior=set_prior("normal(0, 5)", class="b"))## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Start sampling## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: MeanToeArea ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## Data: lizards (Number of observations: 163)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.89 1.65 -4.03 2.39 1.00 2045 2181
## SVL_cm 0.04 0.65 -1.23 1.29 1.00 2109 2157
## HurricaneAfter -0.54 0.30 -1.14 0.05 1.00 2318 2359
## OriginWaterCay -0.08 0.32 -0.71 0.54 1.00 1944 2617
## ISVL_cmE2 0.12 0.06 0.00 0.25 1.00 2185 2152
## SVL_cm:HurricaneAfter 0.16 0.06 0.04 0.28 1.00 2273 2387
## SVL_cm:OriginWaterCay 0.01 0.06 -0.12 0.14 1.00 1926 2531
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.24 0.01 0.21 0.27 1.00 3658 2551
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(conditional_effects(brm_lm4, effects="SVL_cm:Hurricane", conditions=data.frame(Origin=levels(lizards$Origin), cond__=levels(lizards$Origin)), plot=FALSE), facet_args=list(ncol=2), points=TRUE)
Back to the questions
Goals
- describe how traits differ before/after hurricanes
- account for confounding factor: overall size
- also possibly account for differences between islands and sexes
For instance, we’d like to say something like:
Remaining male lizards of average size on Pine Cay had toepads that were X% larger by area (95% CI of L-U%). This represents roughly Y standard deviations of toepad area. Toepad area did not differ substantially by island, or by sex after accounting for length.
But, we have twelve variables! What about the rest?
We could do twelve univariate analyses.
But, what if (say) toe length increases only because it’s correlated with toepad area?
Multivariate
Here’s the same model, but with all twelve responses.
brm_mvlm <- brm(
brmsformula(mvbind(Femur, Tibia, Metatarsal, LongestToe, Humerus,
Radius, Metacarpal, LongestFinger, FingerCount,
ToeCount, MeanFingerArea, MeanToeArea)
~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)) + set_rescor(TRUE),
data=lizards, prior=set_prior("normal(0, 5)", class="b"))## Warning: Rows containing NAs were excluded from the model.## Warning: Specifying global priors for regression coefficients in multivariate models is deprecated. Please specify priors separately for each response variable.## Compiling Stan program...## Start sampling## Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problems## Family: MV(gaussian, gaussian, gaussian, gaussian, gaussian, gaussian, gaussian, gaussian, gaussian, gaussian, gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: Femur ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## Tibia ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## Metatarsal ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## LongestToe ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## Humerus ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## Radius ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## Metacarpal ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## LongestFinger ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## FingerCount ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## ToeCount ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## MeanFingerArea ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## MeanToeArea ~ SVL_cm * (Hurricane + Origin) + I(SVL_cm^2)
## Data: lizards (Number of observations: 163)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Femur_Intercept -0.61 3.96 -8.03 7.35 1.00 1326 1942
## Tibia_Intercept -3.64 2.94 -9.40 1.90 1.00 1141 2003
## Metatarsal_Intercept -4.52 2.01 -8.54 -0.73 1.01 1159 2115
## LongestToe_Intercept -2.62 2.84 -8.26 2.95 1.00 1403 2337
## Humerus_Intercept -1.51 3.11 -7.69 4.59 1.00 1308 2422
## Radius_Intercept -0.75 1.74 -4.06 2.74 1.01 1057 1442
## Metacarpal_Intercept 1.03 1.73 -2.36 4.40 1.00 1962 2531
## LongestFinger_Intercept 1.32 2.09 -2.62 5.54 1.00 1392 2016
## FingerCount_Intercept -1.67 5.12 -11.61 8.46 1.00 1939 2670
## ToeCount_Intercept 2.72 5.68 -8.17 13.86 1.00 1799 2621
## MeanFingerArea_Intercept 0.43 1.19 -1.84 2.79 1.00 1560 2264
## MeanToeArea_Intercept -0.51 1.66 -3.76 2.75 1.00 1387 2183
## Femur_SVL_cm 2.56 1.56 -0.56 5.47 1.00 1417 2312
## Femur_HurricaneAfter -1.15 0.76 -2.63 0.38 1.00 1290 1807
## Femur_OriginWaterCay -2.33 0.83 -3.99 -0.71 1.00 1327 2179
## Femur_ISVL_cmE2 -0.03 0.15 -0.32 0.27 1.00 1506 2579
## Femur_SVL_cm:HurricaneAfter 0.13 0.15 -0.18 0.43 1.00 1326 1826
## Femur_SVL_cm:OriginWaterCay 0.45 0.17 0.13 0.79 1.00 1341 2045
## Tibia_SVL_cm 3.82 1.16 1.65 6.07 1.00 1189 2018
## Tibia_HurricaneAfter 0.06 0.55 -1.02 1.12 1.00 1137 1994
## Tibia_OriginWaterCay -0.65 0.61 -1.82 0.54 1.01 1052 1813
## Tibia_ISVL_cmE2 -0.13 0.11 -0.35 0.08 1.00 1243 1995
## Tibia_SVL_cm:HurricaneAfter -0.03 0.11 -0.24 0.19 1.00 1154 1909
## Tibia_SVL_cm:OriginWaterCay 0.15 0.12 -0.09 0.39 1.01 1047 1609
## Metatarsal_SVL_cm 3.39 0.79 1.90 4.96 1.01 1210 2085
## Metatarsal_HurricaneAfter 0.16 0.38 -0.60 0.90 1.00 943 2116
## Metatarsal_OriginWaterCay -0.08 0.42 -0.89 0.77 1.01 1099 1871
## Metatarsal_ISVL_cmE2 -0.19 0.08 -0.34 -0.04 1.01 1264 2066
## Metatarsal_SVL_cm:HurricaneAfter -0.04 0.08 -0.19 0.11 1.00 949 2168
## Metatarsal_SVL_cm:OriginWaterCay 0.00 0.08 -0.17 0.17 1.01 1104 1889
## LongestToe_SVL_cm 2.47 1.12 0.26 4.70 1.00 1447 2572
## LongestToe_HurricaneAfter 0.48 0.56 -0.64 1.58 1.00 1303 2170
## LongestToe_OriginWaterCay 0.17 0.59 -0.97 1.32 1.00 1383 2101
## LongestToe_ISVL_cmE2 -0.09 0.11 -0.30 0.13 1.00 1501 2551
## LongestToe_SVL_cm:HurricaneAfter -0.15 0.11 -0.37 0.08 1.00 1321 2208
## LongestToe_SVL_cm:OriginWaterCay -0.06 0.12 -0.29 0.17 1.00 1391 1967
## Humerus_SVL_cm 2.19 1.22 -0.22 4.64 1.00 1347 2339
## Humerus_HurricaneAfter -0.77 0.59 -1.97 0.39 1.00 1178 1914
## Humerus_OriginWaterCay -0.13 0.66 -1.40 1.18 1.00 1432 2273
## Humerus_ISVL_cmE2 -0.02 0.12 -0.26 0.21 1.00 1399 2490
## Humerus_SVL_cm:HurricaneAfter 0.22 0.12 -0.02 0.46 1.00 1176 2092
## Humerus_SVL_cm:OriginWaterCay 0.05 0.13 -0.22 0.30 1.00 1432 2217
## Radius_SVL_cm 1.94 0.68 0.57 3.23 1.01 1107 1516
## Radius_HurricaneAfter -0.65 0.33 -1.29 -0.01 1.00 1056 2080
## Radius_OriginWaterCay -0.55 0.37 -1.27 0.16 1.01 1065 1894
## Radius_ISVL_cmE2 -0.04 0.07 -0.17 0.09 1.01 1168 1777
## Radius_SVL_cm:HurricaneAfter 0.13 0.07 0.01 0.26 1.00 1039 2121
## Radius_SVL_cm:OriginWaterCay 0.13 0.07 -0.02 0.27 1.01 1063 1936
## Metacarpal_SVL_cm 0.22 0.68 -1.12 1.55 1.00 2036 2500
## Metacarpal_HurricaneAfter -0.02 0.32 -0.65 0.63 1.00 2294 2928
## Metacarpal_OriginWaterCay -0.30 0.35 -0.98 0.38 1.00 2132 2835
## Metacarpal_ISVL_cmE2 0.02 0.07 -0.11 0.15 1.00 2122 2522
## Metacarpal_SVL_cm:HurricaneAfter 0.02 0.06 -0.11 0.14 1.00 2265 3011
## Metacarpal_SVL_cm:OriginWaterCay 0.09 0.07 -0.05 0.22 1.00 2142 2845
## LongestFinger_SVL_cm 0.31 0.82 -1.34 1.85 1.00 1433 2128
## LongestFinger_HurricaneAfter -0.40 0.39 -1.18 0.35 1.00 1388 2206
## LongestFinger_OriginWaterCay -0.88 0.42 -1.70 -0.05 1.00 1441 2178
## LongestFinger_ISVL_cmE2 0.04 0.08 -0.12 0.20 1.00 1491 2145
## LongestFinger_SVL_cm:HurricaneAfter 0.09 0.08 -0.06 0.25 1.00 1361 2080
## LongestFinger_SVL_cm:OriginWaterCay 0.19 0.08 0.02 0.35 1.00 1429 2370
## FingerCount_SVL_cm 4.17 2.02 0.20 8.03 1.00 2025 2674
## FingerCount_HurricaneAfter 1.27 1.08 -0.85 3.32 1.00 1160 2470
## FingerCount_OriginWaterCay -0.43 1.09 -2.56 1.74 1.00 1741 2517
## FingerCount_ISVL_cmE2 -0.31 0.20 -0.68 0.09 1.00 2117 2790
## FingerCount_SVL_cm:HurricaneAfter -0.17 0.22 -0.59 0.25 1.00 1163 2513
## FingerCount_SVL_cm:OriginWaterCay 0.11 0.22 -0.33 0.54 1.00 1763 2489
## ToeCount_SVL_cm 3.31 2.24 -1.09 7.65 1.00 1883 2754
## ToeCount_HurricaneAfter -0.18 1.21 -2.55 2.15 1.00 1740 2588
## ToeCount_OriginWaterCay -1.58 1.25 -4.02 0.87 1.00 1709 2389
## ToeCount_ISVL_cmE2 -0.21 0.22 -0.64 0.21 1.00 1980 2772
## ToeCount_SVL_cm:HurricaneAfter 0.08 0.24 -0.40 0.55 1.00 1730 2445
## ToeCount_SVL_cm:OriginWaterCay 0.34 0.25 -0.16 0.83 1.00 1722 2161
## MeanFingerArea_SVL_cm -0.33 0.47 -1.25 0.56 1.00 1656 2487
## MeanFingerArea_HurricaneAfter -0.30 0.22 -0.74 0.13 1.00 1517 2585
## MeanFingerArea_OriginWaterCay -0.45 0.23 -0.89 0.02 1.00 1599 2407
## MeanFingerArea_ISVL_cmE2 0.11 0.05 0.02 0.20 1.00 1767 2484
## MeanFingerArea_SVL_cm:HurricaneAfter 0.10 0.04 0.02 0.19 1.00 1492 2493
## MeanFingerArea_SVL_cm:OriginWaterCay 0.11 0.05 0.01 0.20 1.00 1599 2448
## MeanToeArea_SVL_cm -0.11 0.65 -1.40 1.16 1.00 1442 2101
## MeanToeArea_HurricaneAfter -0.54 0.30 -1.11 0.05 1.00 1526 2728
## MeanToeArea_OriginWaterCay -0.11 0.33 -0.74 0.54 1.00 1456 1979
## MeanToeArea_ISVL_cmE2 0.14 0.06 0.02 0.26 1.00 1513 2167
## MeanToeArea_SVL_cm:HurricaneAfter 0.16 0.06 0.04 0.28 1.00 1440 2749
## MeanToeArea_SVL_cm:OriginWaterCay 0.02 0.07 -0.11 0.14 1.00 1437 2107
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_Femur 0.62 0.03 0.55 0.69 1.00 3343 2938
## sigma_Tibia 0.45 0.02 0.40 0.50 1.00 3088 3737
## sigma_Metatarsal 0.31 0.02 0.28 0.35 1.00 3250 3683
## sigma_LongestToe 0.44 0.02 0.39 0.49 1.00 4179 3210
## sigma_Humerus 0.48 0.03 0.43 0.54 1.00 3421 3249
## sigma_Radius 0.26 0.01 0.24 0.29 1.00 3533 3181
## sigma_Metacarpal 0.25 0.01 0.23 0.28 1.00 4424 2628
## sigma_LongestFinger 0.31 0.02 0.28 0.35 1.00 4351 3265
## sigma_FingerCount 0.86 0.05 0.77 0.96 1.00 3889 3272
## sigma_ToeCount 0.99 0.06 0.89 1.11 1.00 3635 3197
## sigma_MeanFingerArea 0.17 0.01 0.15 0.19 1.00 4317 2736
## sigma_MeanToeArea 0.24 0.01 0.22 0.27 1.00 4779 2982
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(Femur,Tibia) 0.46 0.06 0.33 0.57 1.00 3071 2958
## rescor(Femur,Metatarsal) 0.29 0.07 0.14 0.42 1.00 2603 2811
## rescor(Tibia,Metatarsal) 0.61 0.05 0.51 0.70 1.00 3177 3052
## rescor(Femur,LongestToe) 0.15 0.07 0.01 0.29 1.00 3759 3338
## rescor(Tibia,LongestToe) 0.20 0.07 0.06 0.34 1.00 3806 2911
## rescor(Metatarsal,LongestToe) 0.39 0.07 0.26 0.52 1.00 4009 3279
## rescor(Femur,Humerus) 0.32 0.07 0.18 0.46 1.00 3337 2803
## rescor(Tibia,Humerus) 0.44 0.06 0.30 0.55 1.00 3686 3274
## rescor(Metatarsal,Humerus) 0.41 0.06 0.28 0.53 1.00 3943 3199
## rescor(LongestToe,Humerus) 0.25 0.07 0.10 0.39 1.00 5297 3405
## rescor(Femur,Radius) 0.47 0.06 0.34 0.58 1.00 3059 2990
## rescor(Tibia,Radius) 0.54 0.06 0.42 0.64 1.00 3551 3069
## rescor(Metatarsal,Radius) 0.50 0.06 0.37 0.60 1.00 3646 3341
## rescor(LongestToe,Radius) 0.27 0.07 0.13 0.40 1.00 4675 3416
## rescor(Humerus,Radius) 0.49 0.06 0.37 0.60 1.00 3712 3051
## rescor(Femur,Metacarpal) 0.04 0.08 -0.12 0.19 1.00 4359 3110
## rescor(Tibia,Metacarpal) 0.12 0.07 -0.03 0.26 1.00 4481 3051
## rescor(Metatarsal,Metacarpal) 0.14 0.07 -0.01 0.28 1.00 5011 3271
## rescor(LongestToe,Metacarpal) 0.24 0.07 0.08 0.37 1.00 5055 2774
## rescor(Humerus,Metacarpal) 0.21 0.07 0.06 0.36 1.00 5162 3194
## rescor(Radius,Metacarpal) 0.15 0.07 0.01 0.30 1.00 5331 3037
## rescor(Femur,LongestFinger) 0.15 0.08 -0.00 0.29 1.00 3492 3222
## rescor(Tibia,LongestFinger) 0.23 0.07 0.09 0.36 1.00 3935 3342
## rescor(Metatarsal,LongestFinger) 0.35 0.07 0.22 0.48 1.00 4523 3392
## rescor(LongestToe,LongestFinger) 0.39 0.07 0.26 0.51 1.00 4945 2974
## rescor(Humerus,LongestFinger) 0.23 0.07 0.09 0.37 1.00 4881 2866
## rescor(Radius,LongestFinger) 0.22 0.07 0.07 0.35 1.00 4460 3452
## rescor(Metacarpal,LongestFinger) 0.20 0.07 0.06 0.35 1.00 4700 3198
## rescor(Femur,FingerCount) 0.19 0.07 0.03 0.33 1.00 3951 3515
## rescor(Tibia,FingerCount) 0.24 0.07 0.09 0.38 1.00 3968 3522
## rescor(Metatarsal,FingerCount) 0.28 0.07 0.14 0.41 1.00 4626 3423
## rescor(LongestToe,FingerCount) 0.19 0.07 0.04 0.33 1.00 4608 3377
## rescor(Humerus,FingerCount) 0.11 0.07 -0.04 0.25 1.00 4790 3193
## rescor(Radius,FingerCount) 0.23 0.07 0.09 0.37 1.00 4479 3513
## rescor(Metacarpal,FingerCount) 0.10 0.08 -0.05 0.25 1.00 4190 3316
## rescor(LongestFinger,FingerCount) 0.24 0.07 0.10 0.38 1.00 4066 3213
## rescor(Femur,ToeCount) 0.11 0.08 -0.04 0.26 1.00 4112 3387
## rescor(Tibia,ToeCount) 0.29 0.07 0.15 0.42 1.00 4253 3041
## rescor(Metatarsal,ToeCount) 0.22 0.07 0.08 0.35 1.00 5126 3383
## rescor(LongestToe,ToeCount) 0.10 0.08 -0.05 0.25 1.00 5303 2992
## rescor(Humerus,ToeCount) 0.12 0.08 -0.03 0.27 1.00 5187 3413
## rescor(Radius,ToeCount) 0.13 0.07 -0.02 0.27 1.00 5190 3620
## rescor(Metacarpal,ToeCount) 0.09 0.08 -0.06 0.24 1.00 4425 3403
## rescor(LongestFinger,ToeCount) 0.22 0.07 0.07 0.36 1.00 4430 2754
## rescor(FingerCount,ToeCount) 0.46 0.06 0.34 0.58 1.00 4383 2974
## rescor(Femur,MeanFingerArea) 0.19 0.07 0.04 0.33 1.00 4947 3133
## rescor(Tibia,MeanFingerArea) 0.15 0.07 -0.00 0.29 1.00 4630 3087
## rescor(Metatarsal,MeanFingerArea) 0.12 0.07 -0.02 0.27 1.00 5123 3311
## rescor(LongestToe,MeanFingerArea) 0.24 0.07 0.09 0.38 1.00 4666 3171
## rescor(Humerus,MeanFingerArea) 0.14 0.07 -0.01 0.27 1.00 4556 3253
## rescor(Radius,MeanFingerArea) 0.18 0.07 0.04 0.32 1.00 5263 3248
## rescor(Metacarpal,MeanFingerArea) 0.15 0.08 -0.01 0.30 1.00 3866 3311
## rescor(LongestFinger,MeanFingerArea) 0.21 0.07 0.07 0.35 1.00 3887 3172
## rescor(FingerCount,MeanFingerArea) 0.31 0.07 0.17 0.45 1.00 4302 3059
## rescor(ToeCount,MeanFingerArea) 0.16 0.08 0.01 0.31 1.00 4659 3182
## rescor(Femur,MeanToeArea) 0.19 0.07 0.05 0.33 1.00 4033 3221
## rescor(Tibia,MeanToeArea) 0.24 0.07 0.10 0.37 1.00 4662 3577
## rescor(Metatarsal,MeanToeArea) 0.20 0.07 0.05 0.34 1.00 4316 3453
## rescor(LongestToe,MeanToeArea) 0.30 0.07 0.17 0.44 1.00 4450 3180
## rescor(Humerus,MeanToeArea) 0.21 0.07 0.06 0.35 1.00 4474 3492
## rescor(Radius,MeanToeArea) 0.36 0.07 0.22 0.48 1.00 4670 3427
## rescor(Metacarpal,MeanToeArea) 0.18 0.07 0.03 0.32 1.00 4181 3032
## rescor(LongestFinger,MeanToeArea) 0.30 0.07 0.16 0.43 1.00 4509 3483
## rescor(FingerCount,MeanToeArea) 0.05 0.08 -0.09 0.20 1.00 4099 3522
## rescor(ToeCount,MeanToeArea) 0.14 0.07 -0.01 0.28 1.00 4510 3349
## rescor(MeanFingerArea,MeanToeArea) 0.45 0.06 0.33 0.57 1.00 4497 2858
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Wait, is “longest toe” different or not?
For concreteness, let’s look at the posterior distribution of SVL=5cm lizards on Pine Cay, both (a) the posterior predictive distribution (reflecting variation in the population), and (b) the posterior distribution of mean measurements (reflecting uncertainty).
pred_df <- expand.grid(Origin="Pine Cay", SVL_cm=5, Hurricane=levels(lizards$Hurricane))
post_pred_array <- posterior_predict(brm_mvlm, newdata=pred_df, resp=responses)
post_mean_array <- posterior_epred(brm_mvlm, newdata=pred_df, resp=responses)
dimnames(post_pred_array)[[2]] <- dimnames(post_mean_array)[[2]] <- levels(lizards$Hurricane)
# distribution of lizard measurements
post_pred <- rbind(cbind(data.frame(post_pred_array[,"Before",]), Hurricane="Before"),
cbind(data.frame(post_pred_array[,"After",]), Hurricane="After"))
# posterior distribution of mean lizard measurements
post_mean <- rbind(cbind(data.frame(post_mean_array[,"Before",]), Hurricane="Before"),
cbind(data.frame(post_mean_array[,"After",]), Hurricane="After"))Posterior predictive distribution of responses: 
Posterior distribution of means: 
Now, we’ll summarize that by getting posterior means and 95% credible intervals for each of the differences, as a percent of the “before” value.
ldf <- na.omit(lizards[,setdiff(colnames(lizards), responses)])
ldf <- rbind(ldf, ldf)
ldf$Hurricane <- rep(levels(lizards$Hurricane), each=nrow(ldf)/2)
pldf <- posterior_epred(brm_mvlm, newdata=ldf, summarise=FALSE)
diff_ldf <- apply((pldf[,ldf$Hurricane == "After",] - pldf[,ldf$Hurricane == "Before",])/pldf[,ldf$Hurricane == "Before",], c(1,3), mean)
stopifnot(all(colnames(diff_ldf) == responses))
post_table <- data.frame("mean difference"=100 * colMeans(diff_ldf),
"lower 95% CI"=100 * colQuantiles(diff_ldf, probs=0.025),
"upper 95% CI"=100 * colQuantiles(diff_ldf, probs=1-0.025),
check.names=FALSE)## mean difference lower 95% CI upper 95% CI
## Femur -4.78348 -6.5020 -3.0562
## Tibia -0.55969 -1.7667 0.6545
## Metatarsal -0.43780 -1.7759 0.9297
## LongestToe -3.32220 -5.1501 -1.4163
## Humerus 3.23052 1.3556 5.0580
## Radius 0.05234 -1.0505 1.1239
## Metacarpal 2.45169 -0.6634 5.4179
## LongestFinger 0.61358 -1.9592 3.3166
## FingerCount 3.77986 1.3584 6.2031
## ToeCount 1.38851 -0.8867 3.7258
## MeanFingerArea 14.22498 9.5351 18.9842
## MeanToeArea 10.46696 6.7681 14.4858
Back to the questions
Goals
- describe how traits differ before/after hurricanes
- account for confounding factor: overall size
- also possibly account for differences between islands and sexes