3
$\begingroup$

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 

enter image description here

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

enter image description here

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

enter image description here

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()
$\endgroup$
5
  • 2
    $\begingroup$ Can you give some details about your Area variable? What does it represent and what values can it assume? In particular, how many unique values does Area have? Also, what fraction of values - if any - of ProportionBirdsScavenging consist of 0 or 1 values? While we're at it, can you say more about OverheadCover? $\endgroup$ – Isabella Ghement Jan 23 '20 at 16:47
  • 2
    $\begingroup$ When you say the reliability of the p-value, which one? You're showing quite a few which aren't comparable (some are tests of fixed effects of single terms, some are special tests for ranefs, some are GLRT comparing models that differ by several terms). Besides providing some info on exactly what you are doing --- just showing model code and plots without the code to generate them is not very useful --- do read ?summary.gam and ?anova.gam as 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:44
  • $\begingroup$ @ Isabella Ghement Please see my edit. I should probably have explained more about my data at first instance, my apologies. @ Reinstate Monica I'm talking about the all pvalues above (summary, joint_tests, and lrtest) as being unreliable/difficult to interpret, because the OverheadCover still 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
  • $\begingroup$ How did you create your random data? $\endgroup$ – Gavin Simpson Jan 23 '20 at 23:38
  • $\begingroup$ I simply used the =RAND() function in excel to create a column with random numbers between 0 and 1 for x (OverheadCover). The y-values (ProportionBirdsScavenging) remain the same. Please see my edit, as I added data for a reconstruction. $\endgroup$ – Peter Jan 24 '20 at 7:46
0
$\begingroup$

What happens is that the regression forces a line through the heavily weighted data points. The solution is to scale the weights to mitigate their effect. Analysing the data with scaled weights yields reliable p-values.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.