I'm investigating the effect of overhead cover (tree canopy) on the proportion of birds scavenging on carcasses left out in nature. These data are taken from different national parks with different habitats and species pool, hence I am including Area as random effect. I'm running a generalised additive model mgcv::gam with betar family. See data for a reconstruction below.
I want to check wether OverheadCover has a significant effect on PropBirdsScavenging. In first instance I took the P-value of the summary function. However, after experimenting with randomly generated data (made using the =RAND() function in excel), I doubt if this is a reliable p-value to check for significance of OverheadCover (see below). After some research online, I found that besides the summary() function, the emmeans ::joint_tests() and lmtest::lrtest() could also be used, so I included them in the examples below in order to provide extra information and clarify the problem. Which of these tests are reliable in my case, if any? All the plots below show the 95% confidence interval.
> mygam <- mgcv::gam(ProportionBirdsScavenging ~ OverheadCover + s(Area, bs="re"), family=betar(link="logit"), data = df_prop_birds_eating, weights = pointWeight)
> summary(mygam)
Family: Beta regression(5.549)
Link function: logit
Formula:
ProportionBirdsScavenging ~ OverheadCover + s(Area, bs = "re")
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.7564 0.1171 15.00 <2e-16 ***
OverheadCover -4.3864 0.1040 -42.16 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(Area) 4.799 5 196 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.777 Deviance explained = 80.6%
-REML = -3628.8 Scale est. = 1 n = 35
> lrtest(mygam)
Likelihood ratio test
Model 1: ProportionBirdsScavenging ~ OverheadCover + s(Area, bs = "re")
Model 2: ProportionBirdsScavenging ~ 1
#Df LogLik Df Chisq Pr(>Chisq)
1 7.9873 3645.6
2 2.0000 959.6 -5.9873 5371.9 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> joint_tests(mygam)
model term df1 df2 F.ratio p.value
Area 5 28.2 39.205 <.0001
So far, so good. Here I wanted to check the reliability of the p-values, so I created some random data and rerun the analysis, expecting to get a very large p-value.
> mygam <- mgcv::gam(ProportionBirdsScavenging ~ OverheadCover + s(Area, bs="re"), family=betar(link="logit"), data = df_prop_birds_eating, weights = pointWeight)
> summary(mygam)
Family: Beta regression(3.61)
Link function: logit
Formula:
ProportionBirdsScavenging ~ OverheadCover + s(Area, bs = "re")
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.38816 0.47704 -0.814 0.416
OverheadCover -0.26252 0.05028 -5.221 1.78e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(Area) 4.99 5 6033 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.617 Deviance explained = 69.3%
-REML = -2841.6 Scale est. = 1 n = 35
> lrtest(mygam)
Likelihood ratio test
Model 1: ProportionBirdsScavenging ~ OverheadCover + s(Area, bs = "re")
Model 2: ProportionBirdsScavenging ~ 1
#Df LogLik Df Chisq Pr(>Chisq)
1 8 2866.06
2 2 959.62 -6 3812.9 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> joint_tests(mygam)
model term df1 df2 F.ratio p.value
Area 5 28.01 1206.622 <.0001
The confidence interval seems OK, but the p-values < 0.001 surprised me, especially when looking at the plot. The effect of X should not be significant here. I thought maybe I'll try it without the random effect.
> mygam <- mgcv::gam(ProportionBirdsScavenging ~ OverheadCover, family=betar(link="logit"), data = df_prop_birds_eating, weights = pointWeight)
> summary(mygam)
Family: Beta regression(1.334)
Link function: logit
Formula:
ProportionBirdsScavenging ~ OverheadCover
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.48564 0.03424 -14.18 <2e-16 ***
OverheadCover -0.01284 0.05842 -0.22 0.826
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = -0.0351 Deviance explained = 0.884%
-REML = -954.63 Scale est. = 1 n = 35
> lrtest(mygam)
Likelihood ratio test
Model 1: ProportionBirdsScavenging ~ OverheadCover
Model 2: ProportionBirdsScavenging ~ 1
#Df LogLik Df Chisq Pr(>Chisq)
1 3 959.64
2 2 959.62 -1 0.0491 0.8246
> joint_tests(mygam)
# only white space
Here we see that the p-value is very insignificant (as I expect it to be), but the confidence interval is very small.
Does somebody know what is going on here? Why are these p-values unreliable with a random effect, and why is the confidence interval unreliable without a random effect? Am I interpreting this wrong, is there maybe a different test to find reliable p-values in my particular case?
Please see the script and data below for a reconstruction. Unfortunately, I could not retrieve the exact same random data as the plots above, but the idea is the same.
library(mgcv)
library(sjstats)
library(emmeans)
library(lmtest)
library(dplyr)
library(ggplot2)
rm(list = ls())
df_prop_birds_eating <- data.frame(ProportionBirdsScavenging = c(0.661904761904762, 0.406192519926425, 0.694285714285714, 0.0142857142857143, 0.0142857142857143, 0.0584415584415584, 0.257142857142857, 0.701154401154401, 0.625521669341894, 0.871428571428571, 0.930204081632653, 0.452995391705069, 0.0483709273182957, 0.5, 0.484821428571429, 0.0142857142857143, 0.0345238095238095, 0.248148148148148, 0.795866569978245, 0.661904761904762, 0.985714285714286, 0.985714285714286, 0.985714285714286, 0.550746268656716, 0.208571428571429, 0.0142857142857143, 0.0142857142857143, 0.0142857142857143, 0.0142857142857143, 0.0142857142857143, 0.14030888030888, 0.101142857142857, 0.573982970671712, 0.0142857142857143, 0.0142857142857143),
pointWeight = c(3, 233, 10, 89, 4, 22, 44, 99, 89, 17, 35, 341, 57, 36, 128, 39, 144, 54, 394, 12, 46, 229, 55, 67, 5, 28, 2, 160, 124, 294, 555, 425, 302, 116, 48),
OverheadCover = c(0.7, 0.671, 0.6795, 0.79, 0.62, 0.62, 0.6413, 0.089, 0.4603, 0.04, 0.0418, 0.46, 0.5995, 0.532, 0.65, 0.6545, 0.74, 0.74, 0.02, 0.02, 0, 0, 0, 0.45, 0.8975, 0.92, 0.898, 0.89, 0.86, 0.69, 0.755, 0.775, 0.585, 0.585, 0.55),
Random_OverheadCover = c(0.413296794, 0.347613463, 0.349795504, 0.961276169, 0.95789013, 0.035865447, 0.324590684, 0.44398489, 0.199835127, 0.53304746, 0.03317392, 0.9666195, 0.306412904, 0.208397613, 0.461614514, 0.209814992, 0.136107136, 0.107251692, 0.742954177, 0.620547648, 0.724865263, 0.290264773, 0.877152686, 0.26345651, 0.113437669, 0.333750829, 0.167373418, 0.431189319, 0.158608563, 0.026153635, 0.443326766, 0.450652285, 0.787105741, 0.833940574, 0.040695708),
Area = c("Markiezaat", "Hamert", "Hamert", "Hamert", "Hamert", "Hamert", "Hamert", "Hamert", "Hamert", "KempenBroek", "KempenBroek", "KempenBroek", "KempenBroek", "KempenBroek", "KempenBroek", "KempenBroek", "KempenBroek", "KempenBroek", "Markiezaat", "Markiezaat", "Markiezaat", "Markiezaat", "Markiezaat", "Meinweg", "Meinweg", "Meinweg", "PlankenWambuis", "PlankenWambuis", "PlankenWambuis", "PlankenWambuis", "PlankenWambuis", "PlankenWambuis", "Valkenhorst", "Valkenhorst", "KempenBroek"))
mygam <- mgcv::gam(ProportionBirdsScavenging ~ OverheadCover + s(Area, bs="re"), family=betar(link="logit"), data = df_prop_birds_eating, weights = pointWeight)
summary(mygam)
lrtest(mygam)
joint_tests(mygam)
new.xgam <- expand.grid(OverheadCover = seq(min(df_prop_birds_eating$OverheadCover), max(df_prop_birds_eating$OverheadCover), length.out = 1000))
new.xgam$Area <- "Hamert" # pad new.xgam with an arbitrary value for variable Area, otherwise predict() won't work -> https://stackoverflow.com/questions/54411851/mgcv-how-to-use-exclude-argument-in-predict-gam
new.ygam <- data.frame(predict.gam(mygam, new.xgam, type = "response", exclude = "s(Area)", se.fit = TRUE)) %>% rename(ProportionBirdsScavenging = fit, SE = se.fit)
addThesegam <- mutate(data.frame(new.xgam, new.ygam),
lwr = ProportionBirdsScavenging - 1.96 * SE,
upr = ProportionBirdsScavenging + 1.96 * SE)
ggplot(df_prop_birds_eating, aes(x = OverheadCover, y = ProportionBirdsScavenging)) +
geom_point(aes(size = pointWeight), pch = 21, fill = I("red"), alpha = 0.3) +
geom_smooth(aes(ymin = lwr, ymax = upr), data = addThesegam, stat = 'identity') +
theme_bw()
rm(mygam, new.xgam, new.ygam, addThesegam)
mygam <- mgcv::gam(ProportionBirdsScavenging ~ Random_OverheadCover + s(Area, bs="re"), family=betar(link="logit"), data = df_prop_birds_eating, weights = pointWeight)
summary(mygam)
lrtest(mygam)
joint_tests(mygam)
new.xgam <- expand.grid(Random_OverheadCover = seq(min(df_prop_birds_eating$Random_OverheadCover), max(df_prop_birds_eating$Random_OverheadCover), length.out = 1000))
new.xgam$Area <- "Hamert" # pad new.xgam with an arbitrary value for variable Area, otherwise predict() won't work -> https://stackoverflow.com/questions/54411851/mgcv-how-to-use-exclude-argument-in-predict-gam
new.ygam <- data.frame(predict.gam(mygam, new.xgam, type = "response", exclude = "s(Area)", se.fit = TRUE)) %>% rename(ProportionBirdsScavenging = fit, SE = se.fit)
addThesegam <- mutate(data.frame(new.xgam, new.ygam),
lwr = ProportionBirdsScavenging - 1.96 * SE,
upr = ProportionBirdsScavenging + 1.96 * SE)
ggplot(df_prop_birds_eating, aes(x = Random_OverheadCover, y = ProportionBirdsScavenging)) +
geom_point(aes(size = pointWeight), pch = 21, fill = I("red"), alpha = 0.3) +
geom_smooth(aes(ymin = lwr, ymax = upr), data = addThesegam, stat = 'identity') +
theme_bw()
rm(mygam, new.xgam, new.ygam, addThesegam)
mygam <- mgcv::gam(ProportionBirdsScavenging ~ Random_OverheadCover, family=betar(link="logit"), data = df_prop_birds_eating, weights = pointWeight)
summary(mygam)
lrtest(mygam)
joint_tests(mygam)
new.xgam <- expand.grid(Random_OverheadCover = seq(min(df_prop_birds_eating$Random_OverheadCover), max(df_prop_birds_eating$Random_OverheadCover), length.out = 1000))
new.xgam$Area <- "Hamert" # pad new.xgam with an arbitrary value for variable Area, otherwise predict() won't work -> https://stackoverflow.com/questions/54411851/mgcv-how-to-use-exclude-argument-in-predict-gam
new.ygam <- data.frame(predict.gam(mygam, new.xgam, type = "response", exclude = "s(Area)", se.fit = TRUE)) %>% rename(ProportionBirdsScavenging = fit, SE = se.fit)
addThesegam <- mutate(data.frame(new.xgam, new.ygam),
lwr = ProportionBirdsScavenging - 1.96 * SE,
upr = ProportionBirdsScavenging + 1.96 * SE)
ggplot(df_prop_birds_eating, aes(x = Random_OverheadCover, y = ProportionBirdsScavenging)) +
geom_point(aes(size = pointWeight), pch = 21, fill = I("red"), alpha = 0.3) +
geom_smooth(aes(ymin = lwr, ymax = upr), data = addThesegam, stat = 'identity') +
theme_bw()
?summary.gamand?anova.gamas getting p-values for ranefs is not easy and doing GLRTs tricky, esp as you you may/may not be accounting for selection of smoothness params. $\endgroup$ – Gavin Simpson Jan 23 '20 at 17:44summary,joint_tests, andlrtest) as being unreliable/difficult to interpret, because theOverheadCoverstill has a very significant effect even with random data. However, I'm probably misusing some of these tests. Ultimately, I'm looking for a test which gives me a (reliable) p-value for my specific dataset (which will be significant for my data, but not for the randomly generated data). $\endgroup$ – Peter Jan 23 '20 at 21:25