0
$\begingroup$

My question is about statistical interpretation of a two-level categorical fixed effect in a GAM model regarding the difference between two groups. Altough I found information about the construction of such models, for example here:

I could not find concrete guidance about the significance values and their meaning in model summary, particularly regarding the effect of the categorical variable.


I have two time series of soil moisture measurements over time, obtained in two locations: sun and shade.

Here is reproducible R code to recreate the data:

dat = structure(list(day = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 12, 13, 
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 28, 29, 30, 
31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 52, 53, 54, 55, 
56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 
72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 
88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 
103, 104, 105, 106, 107, 0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 12, 13, 
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 28, 29, 30, 
31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 52, 53, 54, 55, 
56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 
72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 
88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 
103, 104, 105, 106, 107), location = c("sun", "sun", "sun", "sun", 
"sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", 
"sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", 
"sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", 
"sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", 
"sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", 
"sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", 
"sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", 
"sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", 
"sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", 
"sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", "sun", 
"sun", "sun", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade", "shade", 
"shade", "shade", "shade", "shade", "shade", "shade"), value = c(0.19688255, 
0.20738015, 0.21831515, 0.2179871, 0.2166749, 0.2133944, 0.2133944, 
0.20563055, 0.2201741, 0.2223611, 0.2166749, 0.21569075, 0.2116448, 
0.20738015, 0.2059586, 0.1989602, 0.19557035, 0.1878065, 0.1832138, 
0.17698085, 0.1622186, 0.14680025, 0.12230585, 0.09628055, 0.07747235, 
0.0935468, 0.0834866, 0.07550405, 0.070802, 0.06588125, 0.06216335, 
0.05735195, 0.0541808, 0.0524312, 0.0502442, 0.04838525, 0.04619825, 
0.04576085, 0.04357385, 0.042371, 0.0316547, 0.0296864, 0.03088925, 
0.0296864, 0.0279368, 0.0279368, 0.02717135, 0.02717135, 0.02498435, 
0.02585915, 0.02717135, 0.02498435, 0.0242189, 0.02585915, 0.02585915, 
0.02498435, 0.0242189, 0.0242189, 0.0242189, 0.02498435, 0.0242189, 
0.0220319, 0.0237815, 0.0207197, 0.0207197, 0.0189701, 0.0176579, 
0.0207197, 0.0198449, 0.0185327, 0.0176579, 0.0189701, 0.0220319, 
0.0207197, 0.0189701, 0.0207197, 0.0215945, 0.0185327, 0.0198449, 
0.0189701, 0.0207197, 0.0207197, 0.0220319, 0.0189701, 0.0189701, 
0.0189701, 0.01536155, 0.01623635, 0.01623635, 0.0176579, 0.0185327, 
0.0189701, 0.0185327, 0.0185327, 0.0189701, 0.0185327, 0.2315, 
0.2274, 0.2241, 0.2216, 0.2174, 0.2104, 0.2061, 0.2008, 0.2143, 
0.2104, 0.2035, 0.1963, 0.1898, 0.182, 0.1783, 0.1738, 0.1671, 
0.1615, 0.1514, 0.147, 0.1416, 0.1357, 0.1317, 0.1276, 0.1266, 
0.1128, 0.1036, 0.0963, 0.092, 0.0875, 0.0841, 0.0797, 0.0768, 
0.0752, 0.0732, 0.0715, 0.0695, 0.0691, 0.0671, 0.066, 0.0562, 
0.0544, 0.0555, 0.0544, 0.0528, 0.0528, 0.0521, 0.0521, 0.0501, 
0.0509, 0.0521, 0.0501, 0.0494, 0.0509, 0.0509, 0.0501, 0.0494, 
0.0494, 0.0494, 0.0501, 0.0494, 0.0474, 0.049, 0.0462, 0.0462, 
0.0446, 0.0434, 0.0462, 0.0454, 0.0442, 0.0434, 0.0446, 0.0474, 
0.0462, 0.0446, 0.0462, 0.047, 0.0442, 0.0454, 0.0446, 0.0462, 
0.0462, 0.0474, 0.0446, 0.0446, 0.0446, 0.0413, 0.0421, 0.0421, 
0.0434, 0.0442, 0.0446, 0.0442, 0.0442, 0.0446, 0.0442)), class = "data.frame", row.names = c(1L, 
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 12L, 13L, 14L, 15L, 16L, 17L, 
18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 29L, 30L, 31L, 
32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 53L, 
54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 
67L, 68L, 69L, 70L, 71L, 72L, 73L, 74L, 75L, 76L, 77L, 78L, 79L, 
80L, 81L, 82L, 83L, 84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L, 
93L, 94L, 95L, 96L, 97L, 98L, 99L, 100L, 101L, 102L, 103L, 104L, 
105L, 106L, 107L, 108L, 109L, 110L, 111L, 112L, 113L, 114L, 115L, 
116L, 117L, 120L, 121L, 122L, 123L, 124L, 125L, 126L, 127L, 128L, 
129L, 130L, 131L, 132L, 133L, 134L, 135L, 137L, 138L, 139L, 140L, 
141L, 142L, 143L, 144L, 145L, 146L, 147L, 148L, 149L, 150L, 151L, 
161L, 162L, 163L, 164L, 165L, 166L, 167L, 168L, 169L, 170L, 171L, 
172L, 173L, 174L, 175L, 176L, 177L, 178L, 179L, 180L, 181L, 182L, 
183L, 184L, 185L, 186L, 187L, 188L, 189L, 190L, 191L, 192L, 193L, 
194L, 195L, 196L, 197L, 198L, 199L, 200L, 201L, 202L, 203L, 204L, 
205L, 206L, 207L, 208L, 209L, 210L, 211L, 212L, 213L, 214L, 215L, 
216L))

And here is what the first/last six rows look like:

head(dat)
##   day location     value
## 1   0      sun 0.1968825
## 2   1      sun 0.2073802
## 3   2      sun 0.2183152
## 4   3      sun 0.2179871
## 5   4      sun 0.2166749
## 6   5      sun 0.2133944
tail(dat)
##     day location  value
## 211 102    shade 0.0442
## 212 103    shade 0.0446
## 213 104    shade 0.0442
## 214 105    shade 0.0442
## 215 106    shade 0.0446
## 216 107    shade 0.0442

The following plot visualizes the two time series (“sun” location in red, “shade” location in blue).

plot(dat$day[dat$location == "sun"], dat$value[dat$location == "sun"], type = "b", col = "red", ylim = range(dat$value), xlab = "Day", ylab = "value")
points(dat$day[dat$location == "shade"], dat$value[dat$location == "shade"], type = "b", col = "blue")

plot1

My goal is to statistically evaluate whether the “sun” and “shade” location differ, i.e., was soil moisture in the “sun” lower/higher/same compared to “shade” over the study period.

Since the temporal trend is non-linear and wasn’t sure it’s resembling a particular function, I decided to use a GAM as follows:

library(mgcv)

dat$location = factor(dat$location)
fit = gam(value ~ location + s(day, by = location), data = dat, method = "REML")

Following the suggestions in GAM factor smooth interaction--include main effect smooth?

Here is the summary of the model:

summary(fit)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## value ~ location + s(day, by = location)
## 
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0881198  0.0005032  175.12   <2e-16 ***
## locationsun -0.0166002  0.0007116  -23.33   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                        edf Ref.df    F p-value    
## s(day):locationshade 8.364  8.884 1621  <2e-16 ***
## s(day):locationsun   8.889  8.996 2400  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.995   Deviance explained = 99.5%
## -REML = -684.71  Scale est. = 2.4308e-05  n = 192

and the plot:

plot(fit, shade = TRUE, pages = 1, scale = 0, seWithMean = TRUE)

plot2

Based on the locationsun estimate in the Parametric coefficients section (-0.0166002, p<0.001), I conclude that indeed the two series are different, and, specifically, soil moisture was significantly lower in the sun than in the shade.

My question is whether this interpretation is correct? Will be happy to hear any thoughts.

$\endgroup$
  • 2
    $\begingroup$ No, you can't only consider the parametric effects since you fit separate smoothers. Do you have to produce a p-value? Looking at the plot, I'd say no test required since the curves are obviously different. $\endgroup$ – Roland Aug 3 at 11:58
  • $\begingroup$ @Roland Thanks! I see... Then, does it make sense to fit one smoother for both, plus a fixed additive effect that represents how much, on average, one curve is "shifted" from the other? As for your question: no, I don't have to produce p-values, but I have to objectively conclude whether the curves are different or not in some way. Sure, visually they look different, but I need a quantitative measure to support this claim. $\endgroup$ – Michael Dorman Aug 3 at 12:39
  • 2
    $\begingroup$ My first idea was doing fit0 <- gam(value ~ s(day), data = dat, method = "REML"); anova(fit, fit0, test = "Chisq") but then I read help("anova.gam") and developed some doubts. $\endgroup$ – Roland Aug 3 at 13:46

Your Answer

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

Browse other questions tagged or ask your own question.