3
$\begingroup$

Dear StatExchange community,

I am studying disease progression in plant leaves and I am trying to estimate differences between a wild-type and a mutant plant. To achieve this I am using the proportion of healthy pixels in images of plant leaves. Due the type of data I decided to analyze it using gam modeling and beta regression to determinate if the plant genotype has an impact in the disease progression.

So I created two models, one base model with my response variable (size) as a function of hours post infection (hpi). I considered hpi and leaf_id as random factors.

require(mgcv)    
m0 = gam(size ~ s(hpi, k = 4) + s(hpi, leaf_id, bs = "re"), data = data, family = "betar")
m1 = gam(size ~ s(hpi, by = genotype, k = 4) + s(hpi, leaf_id, bs = "re") + genotype, data = data, family = "betar")

All good until here, except when I tested the goodness of fit I did not get any significant p-values even though AIC value improved a lot.

anova(m0, m1, test="Chisq")
Analysis of Deviance Table

Model 1: size ~ s(hpi, k = 4) + s(hpi, leaf_id, bs = "re")
Model 2: size ~ s(hpi, by = genotype, k = 4) + s(hpi, leaf_id, bs = "re") + 
    genotype
  Resid. Df Resid. Dev    Df Deviance Pr(>Chi)
1    92.456     95.862                        
2    89.541     92.029 2.915   3.8333   0.2678
AIC(m0)
[1] -433.6965
AIC(m1)
[1] -451.1168

Do you think there is really not difference between the genotypes or maybe it is incorrect to use anova when one use beta regression?

Here the data I used to fit the models and the fitted model:

leaf_id hpi genotype    size
1_A_01  0   Wild-type   0.999409681
1_A_01  24  Wild-type   0.965263468
1_A_01  48  Wild-type   0.877621373
1_A_01  72  Wild-type   0.527381631
1_A_01  96  Wild-type   0.499122294
1_A_02  0   Wild-type   0.99834779
1_A_02  24  Wild-type   0.9128
1_A_02  48  Wild-type   0.77067519
1_A_02  72  Wild-type   0.419811321
1_A_02  96  Wild-type   0.359519038
1_A_03  0   Wild-type   0.998700455
1_A_03  24  Wild-type   0.941766419
1_A_03  48  Wild-type   0.452631579
1_A_03  72  Wild-type   0.200547778
1_A_03  96  Wild-type   0.119721038
1_A_04  0   Wild-type   0.998331745
1_A_04  24  Wild-type   0.950153555
1_A_04  48  Wild-type   0.7578354
1_A_04  72  Wild-type   0.387827632
1_A_04  96  Wild-type   0.305003427
1_A_05  0   Wild-type   0.994119377
1_A_05  24  Wild-type   0.963669391
1_A_05  48  Wild-type   0.851896813
1_A_05  72  Wild-type   0.618860511
1_A_05  96  Wild-type   0.557251908
1_A_06  0   Wild-type   0.994492044
1_A_06  24  Wild-type   0.93710493
1_A_06  48  Wild-type   0.689569019
1_A_06  72  Wild-type   0.384522734
1_A_06  96  Wild-type   0.309379968
1_A_07  0   Wild-type   0.989961925
1_A_07  24  Wild-type   0.949456174
1_A_07  48  Wild-type   0.809640774
1_A_07  72  Wild-type   0.475455387
1_A_07  96  Wild-type   0.456692913
1_A_08  0   Wild-type   0.979972845
1_A_08  24  Wild-type   0.957887917
1_A_08  48  Wild-type   0.684027778
1_A_08  72  Wild-type   0.472561932
1_A_08  96  Wild-type   0.438386385
1_A_09  0   Wild-type   0.971879106
1_A_09  24  Wild-type   0.956667526
1_A_09  48  Wild-type   0.487238387
1_A_09  72  Wild-type   0.231512071
1_A_09  96  Wild-type   0.179140127
1_A_10  0   Wild-type   0.962017804
1_A_10  24  Wild-type   0.950415592
1_A_10  48  Wild-type   0.607827039
1_A_10  72  Wild-type   0.46189693
1_A_10  96  Wild-type   0.443232662
1_A_11  0   Wild-type   0.968882603
1_A_11  24  Wild-type   0.959830867
1_A_11  48  Wild-type   0.952662722
1_A_11  72  Wild-type   0.88
1_A_11  96  Wild-type   0.764691358
1_A_12  0   Wild-type   0.965517241
1_A_12  24  Wild-type   0.880882353
1_A_12  48  Wild-type   0.313450857
1_A_12  72  Wild-type   0.190072924
1_A_12  96  Wild-type   0.162589928
1_B_01  0   mutant  0.983951856
1_B_01  24  mutant  0.924170616
1_B_01  48  mutant  0.973675964
1_B_01  72  mutant  0.974430068
1_B_01  96  mutant  0.927825261
1_B_02  0   mutant  0.972856418
1_B_02  24  mutant  0.933898305
1_B_02  48  mutant  0.944825953
1_B_02  72  mutant  0.95001179
1_B_02  96  mutant  0.919257652
1_B_03  0   mutant  0.927007299
1_B_03  24  mutant  0.921632653
1_B_03  48  mutant  0.953442879
1_B_03  72  mutant  0.951305025
1_B_03  96  mutant  0.920571882
1_B_04  0   mutant  0.954979129
1_B_04  24  mutant  0.936366254
1_B_04  48  mutant  0.888105727
1_B_04  72  mutant  0.810614525
1_B_04  96  mutant  0.770698849
1_B_05  0   mutant  0.941514075
1_B_05  24  mutant  0.854272864
1_B_05  48  mutant  0.45239486
1_B_05  72  mutant  0.303134886
1_B_05  96  mutant  0.25814978
1_B_06  0   mutant  0.964714389
1_B_06  24  mutant  0.971969069
1_B_06  48  mutant  0.962336245
1_B_06  72  mutant  0.967693971
1_B_06  96  mutant  0.965025374
1_B_07  0   mutant  0.984585742
1_B_07  24  mutant  0.970041841
1_B_07  48  mutant  0.976478275
1_B_07  72  mutant  0.977969974
1_B_07  96  mutant  0.960810143
1_B_08  0   mutant  0.993533584
1_B_08  24  mutant  0.957139955
1_B_08  48  mutant  0.938698833
1_B_08  72  mutant  0.873435055
1_B_08  96  mutant  0.838285145
1_B_09  0   mutant  0.974298565
1_B_09  24  mutant  0.956712673
1_B_09  48  mutant  0.967297408
1_B_09  72  mutant  0.965707287
1_B_09  96  mutant  0.918210198
1_B_10  0   mutant  0.997560976
1_B_10  24  mutant  0.983135392
1_B_10  48  mutant  0.977892756
1_B_10  72  mutant  0.948681397
1_B_10  96  mutant  0.932846715
1_B_11  0   mutant  0.965979637
1_B_11  24  mutant  0.936755047
1_B_11  48  mutant  0.723308979
1_B_11  72  mutant  0.334980005
1_B_11  96  mutant  0.245347698
1_B_12  0   mutant  0.989107764
1_B_12  24  mutant  0.973660714
1_B_12  48  mutant  0.990679095
1_B_12  72  mutant  0.988659106
1_B_12  96  mutant  0.978579481

Fitted model

$\endgroup$
1
$\begingroup$

I wouldn't trust the p-values from the multi-model form of anova.gam(). Reading ?anova.gam Simon gives plenty of warnings about using this, although in the direction opposite to what you see here.

Instead I'd use AIC() as you did already. I would also perhaps compute the differences between the two smooths, and perhaps do a different decomposition but with an ordered genotype:

## pseudo code
ogenotype <- as.ordered(genotype)
m1 <- gam(size ~ s(hpi) + s(hpi, by = ogenotype, k = 4) +
            s(hpi, leaf_id, bs = "re") + ogenotype,
          data = data, family = "betar")

This decomposition will fit a smooth (s(hpi)) for the reference level and a "difference" term relative to the reference level for each other level. This is like the default contrasts in lm() where the intercept is the reference level and the model contains other terms that are the difference between the levels and the reference level.

Then summary(m1) will give you a test directly for a difference between the two.

$\endgroup$
  • $\begingroup$ Thanks for you answer Gavin, indeed the help file suggest not to use p-values for gams model selection. I will stick to AIC then, I tried the second approach you suggested, but apparently I ran out of degrees of freedom. I will try to simulate some data with a higher n. I will declare your answer as definitive. Many thanks for you advice :) $\endgroup$ – Mirko Pavicic Feb 11 at 22:17

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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