3
$\begingroup$

Apologies if a similar question has been asked before.

I have a dataset in which I measured water repellency as a function of water content in soil. I am interested in modeling this relation, as well as the interaction with the addition of a pollutant (in this case copper). Soil samples were incubated for two months at a constant water holding capacity (ranging from 90-20%). In half of these samples, cupper was added. Then water repellency was measured with the molarity of ethanol droplet test (MED). So my response variable is the concentration of ethanol at which a droplet infiltrated a soil surface in a given sample. Generally speaking, the higher the ethanol concentration (in % volume units), the higher the severity of water repellency is.

I am trying to model this data using generalized additive models, because it appears to me that the trend described by the data is non-linear. So I am using GAMs as implemented in the package mgcv in R.

I am struggling to figure out which is the right distribution to assume in this particular case. Since my response variable rep_cat is effectively the percentage of ethanol in water, and the response variable contains lots of 0s.

UPDATE 19.11.22

First I transform my response variable so that it fullfills the conditions of the distributions I have in mind:

rs.cu = rs.cu %>% mutate(rep_cat_b = rep_cat/100) %>% # divides each value by 100 to bound it between 0-1
  mutate(rep_cat_c = rep_cat+0.001) # adds a very small constant to every value, because R considers 0 to be non-positive

Then, I model assuming a binomial distribution and a logit link function:

m1 <- gam(rep_cat_b ~ treatment + s(drought, by=treatment), family = binomial(link = "logit"), 
      data = rs.cu, method = "ML")

The output of the GAM hints that the trend is in fact linear

> summary(m1)

Family: binomial 
Link function: logit 

Formula:
rep_cat_b ~ treatment + s(drought, by = treatment)

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -7.2964     6.9565  -1.049    0.294
treatmentcu  -0.4094    10.1391  -0.040    0.968

Approximate significance of smooth terms:
                        edf Ref.df Chi.sq p-value
s(drought):treatmentcon   1      1  0.230   0.631
s(drought):treatmentcu    1      1  0.346   0.556

R-sq.(adj) =  0.479   Deviance explained = 63.6%
-ML = 1.0518  Scale est. = 1         n = 90

Then I try a second model, but this time I assume a beta distribution:

m2 <- gam(rep_cat_b ~ treatment + s(drought, by=treatment), family = betar(link = "logit"), 
      data = rs.cu, method = "ML")

In this case I get a warning message:

Warning message:
In family$saturated.ll(y, prior.weights, theta) :
saturated likelihood may be inaccurate

And now the results change. The smooth function also indicates the trends are linear. But now the interaction trend is significant:

> summary(m2)

Family: Beta regression(6.771) 
Link function: logit 

Formula:
rep_cat_b ~ treatment + s(drought, by = treatment)

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -5.11580    0.14976 -34.160   <2e-16 ***
treatmentcu  0.04126    0.21181   0.195    0.846    
---
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(drought):treatmentcon   1  1.000  2.470  0.1161  
s(drought):treatmentcu    1  1.001  4.282  0.0386 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.0798   Deviance explained = -0.27%
-ML = -2080.8  Scale est. = 1         n = 90

I fitted both models with method="ML" because I wanted to compare them with AIC, which I did:

> AIC(m1, m2)
         df          AIC
m1 4.000000     9.220146
m2 5.001059 -4151.663549

These results appear to indicate that the second model with beta distribution is a better fit. But still, I am worried about the warning message in the second model and the fact that I have lots of 0s in my response variable.

So I refit the gam using a Gamma distribution with an "identity" link:

m3 <- gam(rep_cat_c ~ treatment + s(drought, by=treatment), family = Gamma(link = "identity"), 
          data = rs.cu, method = "ML", sp=c(0.01, 0.001))

This time, I also get warning messages:

Warning messages:
1: In log(ifelse(y == 0, 1, y/mu)) : NaNs produced
2: Step size truncated due to divergence
3: In log(ifelse(y == 0, 1, y/mu)) : NaNs produced
4: Step size truncated due to divergence
5: In log(ifelse(y == 0, 1, y/mu)) : NaNs produced
6: Step size truncated due to divergence
7: In log(ifelse(y == 0, 1, y/mu)) : NaNs produced
8: Step size truncated due to divergence
9: In log(ifelse(y == 0, 1, y/mu)) : NaNs produced
10: Step size truncated due to divergence
11: In log(ifelse(y == 0, 1, y/mu)) : NaNs produced
12: Step size truncated due to divergence
13: In log(ifelse(y == 0, 1, y/mu)) : NaNs produced
14: Step size truncated due to divergence
15: In log(ifelse(y == 0, 1, y/mu)) : NaNs produced
16: Step size truncated due to divergence

And the model output now suggests non-linear but statistically not significant trends:

Family: Gamma 
Link function: identity 

Formula:
rep_cat_c ~ treatment + s(drought, by = treatment)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   0.4942     0.2770   1.784   0.0787 .
treatmentcu   1.1149     0.8339   1.337   0.1855  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                          edf Ref.df     F p-value
s(drought):treatmentcon 8.522  8.848 0.655   0.745
s(drought):treatmentcu  8.866  8.982 0.468   0.891

R-sq.(adj) =  -1.01   Deviance explained = 76.2%
-ML = -197.28  Scale est. = 3.6284    n = 90

I appreciate any feedback.

Thanks

$\endgroup$
2
  • 1
    $\begingroup$ You mentioned 'error distribution' and 'residual distribution': that's not what you are selecting in the linear model. Rather, you assume a distribution for the conditional expectation. $\endgroup$
    – Firebug
    Nov 11, 2022 at 20:35
  • $\begingroup$ Thanks, @Firebug for pointing this out. I eliminated the words 'residual' or 'error' from the post to avoid misunderstandings. I hope it is clear now what I mean now. $\endgroup$ Nov 12, 2022 at 16:14

2 Answers 2

2
$\begingroup$

From ?formula you'll see that /, in common with other arithmetic operators (+, * etc), has special interpretation. I doubt your (rep_cat / 100) is actually doing what you expect. A negative deviance explained is a sure sign that the fitted model is fundamentally wrong: it suspect that betar() reset all data >= 1 to 1 - eps, which would be any percentage value above 1%.

Refit after converting rep_cat to a proportion outside of the call to gam(), in both cases.

$\endgroup$
1
  • 1
    $\begingroup$ Thanks @Gavin Simpson. I did what you suggested. I transformed my response variable before fitting both models. The output remained identical to the one reported in the original post. $\endgroup$ Nov 12, 2022 at 16:33
2
$\begingroup$

Some thoughts... Looking at the data:

csv <- read.table('https://raw.githubusercontent.com/juan-duenas/datasets/main/rs.cu.csv', sep=',', header=TRUE)
barplot(table(csv$rep_cat))

enter image description here

It seems you have only three values for rep_cat and, even if this is technically percent data, they are far from the upper bound of 100%.

Is this because this is just an example you posted or it is actually the way the data looks like? If the latter, I would consider something like gamma distribution. Also (NB: this is usually a no-no solution): consider binning rep_cat in a small number of factor levels, like "0", "4", "8+" and analyse it using ordinal regression. Finally, I wonder if gam is preferable over glm. I'm not too familiar with gam but it may be that the extra flexibility it provides is not worth it the risk of overfitting.

Maybe I would think about the nature of the data and decide what is a sensible model that generated it rather than going straight to a solution based on diagnostic metrics (I know, not the most useful suggestion.)

$\endgroup$
1
  • $\begingroup$ Thanks '@dariober', these are good points. Please see the updates to the post. Indeed my values are far from the upper bound of 100%. The data is like this, it is not an example. I got a lot of 0s on the higher soil water content categories (very wet to wet soils) because in this condition hydrophobicity is unlikely to appear. I used gams rather than glms because I thought that the trends rather look non-linear. But I accept it might be me forcing my view on the data and thus leading to overfitting. $\endgroup$ Nov 19, 2022 at 14:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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