# Interpreting significance of difference between two GAM smoothed series?

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",
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


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")


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
##
## 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)


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.

• 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. – Roland Aug 3 at 11:58
• @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. – Michael Dorman Aug 3 at 12:39
• 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. – Roland Aug 3 at 13:46