library(tidyr)
library(dplyr)
library(mgcv)
library(gamm4)
library(performance)
library(ggplot2)
# Generate the data
leaf_rand <- data.frame(
Treatment = rep(factor(c("B", "C", "A", "C", "B", "A", "C", "A", "B")),
times = length(2015:2022)),
Year = rep(2015:2022, each = 9),
Plot = factor(rep(1:9, times = 8)),
N.P = round(c(
15.2, 13.6, 15.8, 14.5, 16.1, 17.5, 19.4, 13.8, 14.5,
16.4, 15.7, 14.1, 14.1, 16.1, 18.4, 19.3, 14.8, 16.9,
17.0, 15.1, 18.9, 17.0, 17.3, 17.0, 19.5, 17.0, 16.0,
18.2, 13.1, 19.7, 17.3, 17.3, 17.7, 22.0, 18.3, 14.9,
17.9, 14.4, 16.3, 15.7, 15.8, 20.0, 14.6, 15.4, 12.8,
17.2, 12.4, 13.9, 11.2, 16.3, 22.0, 18.5, 12.4, 12.1,
16.6, 13.6, 17.0, 10.5, 16.2, 18.3, 17.9, 14.5, 14.0,
22.0, 15.4, 20.0, 17.5, 20.8, 21.8, 21.4, 19.4, 17.1
))
)
### GAMM Models ---------------------------------------------------------------
gam_model <- gamm4(N.P ~ Treatment +
s(Year, k = 5, bs = "tp") +
s(Year, Treatment, k = 5, bs = "fs"),
random = ~(1|Plot),
data = leaf_rand, REML = TRUE)
gam_model1 <- gamm(N.P ~ Treatment +
s(Year, k = 5, bs = "tp") +
s(Year, Treatment, k = 5, bs = "fs"),
random = list(Plot=~1),
data = leaf_rand, REML = TRUE)
gam_model2 <- gam(N.P ~ Treatment +
s(Year, k = 5, bs = "tp") +
s(Year, Treatment, k = 5, bs = "fs") +
s(Plot, k = 9, bs = 're'), data = leaf_rand, method="REML")
# Model summaries
summary(gam_model$gam)
summary(gam_model1$gam)
summary(gam_model2)
# Compare the models
compare <- compare_performance(gam_model$gam, gam_model1$gam, gam_model2)
print(compare)
### ggplot by treatment
ggplot(leaf_rand, aes(x = Year, y = N.P, color = Treatment)) +
geom_point() +
geom_smooth(aes(group = Treatment), method = "loess", se = FALSE) +
theme_minimal() +
theme(legend.position = "bottom") +
labs(
x = "Year",
y = "N.P",
color = "Treatment"
)
TheEDIT
Reproducible example with modified models:
gam_model <- gamm4library(N.P ~ stidyr)
library(as.numericdplyr)
library(Yearmgcv), bs = "tp", k = 8
library(gamm4) + Treatment,
random = ~library(1|Plotperformance), data = leaf,
library(ggplot2)
# REMLGenerate =the TRUE)
data
gam_model1leaf_rand <- gamm(Ndata.Pframe(
~ sTreatment = rep(as.numericfactor(Year)c("B", bs"C", ="A", "tp""C", k"B", ="A", 8)"C", +"A", Treatment"B")),
randomtimes = listlength(Plot=~12015:2022)),
data Year = leafrep(2015:2022, methodeach = "REML"9)
,
gam_model2 <- gam(N.PPlot ~= sfactor(as.numericrep(Year), bs = "tp"1:9, ktimes = 8) + Treatment
+s(Plot, bs="re", k=9), data
= leaf, methodN.P = "REML")
The outputs:
> summaryround(gam_model$gam)
c(
Family: gaussian
Link function: identity
Formula:
N15.P ~2, s(as13.numeric(Year)6, bs15.8, =14.5, "tp"16.1, k17.5, =19.4, 13.8) + Treatment
Parametric, coefficients:14.5,
Estimate Std16.4, Error15.7, t14.1, value14.1, Pr(>|t|)16.1, 18.4, 19.3, 14.8, 16.9,
(Intercept) 17.0, 15.2879 1, 18.9, 017.72230, 17.3, 2117.1670, 19.5, <17.0, 2e-16 ***.0,
TreatmentBCN -018.8325 2, 013.4583 -1.817, 19.7, 017.073893, 17.3, 17.7,
TreatmentControl22.0, 18.3, -114.30399,
17.9, 014.45504, 16.3, -215.8667, 15.8, 020.005610, **14.6,
---
Signif15.4, codes:12.8,
0 ‘***’ 017.001 ‘**’2, 012.014, ‘*’13.9, 011.052, ‘16.’3, 022.10, ‘18.5, ’12.4, 12.1
Approximate significance of smooth terms:,
edf Ref16.df F p-value 6,
s(as13.numeric(Year))6, 417.238 0, 410.2385, 316.321 2, 018.01333, *
---
Signif17. codes:9, 14.5, 14.0,
‘***’ 0.001 ‘**’ 022.01 ‘*’ 0.05, ‘15.’4, 020.1 ‘ ’0, 1
R-sq17.(adj) = 5, 020.255 8,
lmer21.REML =8, 29321.79 Scale4, est19. =4, 217.4311 n1
= 72
))
> summary(gam_model1$gam)
Family: gaussian
Link### function:GAMM identityModels
---------------------------------------------------------------
Formula:
N.Pgam_model ~<- sgamm4(asN.numeric(Year), bs = "tp", k = 8)P +~ Treatment
Parametric coefficients:+
Estimate Std. Error t value Pr(>|t|)
s(Intercept)Year, k = 5, bs = "tp") 15.2879+
1.1649 13.123 <2e-16 ***
TreatmentBCN -0.8325 s(Year, Treatment, 1.6475k = -0.5055, bs = "fs"),
0.615
TreatmentControl -1.3039 1.6475 -0.791 0.432 random =
---~(1|Plot),
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
data = leaf_rand, REML = TRUE)
gam_model1 <- gamm(N.P ~ Treatment +
edf Ref.df F p-value
s(as.numeric(Year)) 5.017 , 5.017k 10.39= 5.23e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’, 1
R-sq.(adj)bs = 0.246 "tp") +
Scale est. = 2.4311 n = 72
> summary(gam_model2)
Family: gaussian
Link function: identity
Formula:
N.P ~ s(as.numeric(Year), bs = "tp"Treatment, k = 8) + Treatment + s(Plot5,
bs = "re""fs"),
k = 9)
Parametric coefficients:
random = Estimatelist(Plot=~1),
Std. Error t value Pr(>|t|)
(Intercept) 15.2879 data 1.1649= leaf_rand, 13.123REML = TRUE)
gam_model2 <2e<-16 ***
TreatmentBCNgam(N.P ~ Treatment +
-0.8325 1.6475 -0.505 0.615
TreatmentControl -1.3039s(Year, k = 5, bs 1.6475= "tp") -0.791+
0.432
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘s(Year, ’Treatment, 1
Approximatek significance= of5, smoothbs terms:
= "fs") +
edfs(Plot, Ref.dfk = 9, bs = 're'), Fdata p-value= leaf_rand, method="REML")
# Model summaries
s(as.numericsummary(Yeargam_model$gam)
summary(gam_model1$gam)
summary(gam_model2) 5.017
# 5.927Compare the 8.832models
compare 1.3e<-06 ***compare_performance(gam_model$gam, gam_model1$gam, gam_model2)
sprint(Plotcompare)
### ggplot by treatment
ggplot(leaf_rand, aes(x = Year, y = 5N.552P, color 6.000= 12.397Treatment)) <+
2e-16 ***
---geom_point() +
Signif. codes: geom_smooth(aes(group 0= ‘***’Treatment), 0.001method ‘**’= 0.01"loess", ‘*’se 0.05= ‘.’FALSE) 0.1+
‘ ’theme_minimal() 1
+
R-sq. theme(adj)legend.position = "bottom") 0.651+
labs(
Deviance explained = 71.3%
-REMLx = "Year",
146.9 Scale est.y = 2"N.4311P",
ncolor = 72"Treatment"
)
What is wrong here? Which
Which one should I proceed with? Is
Is it correct to assume that the differences between the treatments are significant according to the $gam parametric coefficient p-values?
It is a first attempt at additive models trying to catch information from an inter annual fluctuating curve, where a linear model perhaps not completely satisfy.
Observed:
Fitted with the first two models:
The overall goal is still to test for differences between the treatments. Along the entire treatment period, as well as within the years.
For the general model structure, this is the initial approach. Looking at things such as adding by=Treatment or a Year interaction, switching the random Plot factor to same penalty structure, etc.
gam_model3 <- gam(N.P ~ s(as.numeric(Year),bs = "tp", k = 8) + as.numeric(Year)*Treatment +
s(as.numeric(Year), Plot, k=8, bs="fs"), data = leaf, method="REML")
Would be happy for any suggestions from the experienced.
Also looking for the appropriate post-hocs.
The models:
gam_model <- gamm4(N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment,
random = ~(1|Plot), data = leaf, REML = TRUE)
gam_model1 <- gamm(N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment,
random = list(Plot=~1), data = leaf, method = "REML")
gam_model2 <- gam(N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment
+s(Plot, bs="re", k=9), data = leaf, method = "REML")
The outputs:
> summary(gam_model$gam)
Family: gaussian
Link function: identity
Formula:
N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.2879 0.7223 21.167 < 2e-16 ***
TreatmentBCN -0.8325 0.4583 -1.817 0.07389 .
TreatmentControl -1.3039 0.4550 -2.866 0.00561 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(as.numeric(Year)) 4.238 4.238 3.321 0.0133 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.255
lmer.REML = 293.79 Scale est. = 2.4311 n = 72
> summary(gam_model1$gam)
Family: gaussian
Link function: identity
Formula:
N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.2879 1.1649 13.123 <2e-16 ***
TreatmentBCN -0.8325 1.6475 -0.505 0.615
TreatmentControl -1.3039 1.6475 -0.791 0.432
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(as.numeric(Year)) 5.017 5.017 10.39 5.23e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.246
Scale est. = 2.4311 n = 72
> summary(gam_model2)
Family: gaussian
Link function: identity
Formula:
N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment + s(Plot,
bs = "re", k = 9)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.2879 1.1649 13.123 <2e-16 ***
TreatmentBCN -0.8325 1.6475 -0.505 0.615
TreatmentControl -1.3039 1.6475 -0.791 0.432
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(as.numeric(Year)) 5.017 5.927 8.832 1.3e-06 ***
s(Plot) 5.552 6.000 12.397 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.651 Deviance explained = 71.3%
-REML = 146.9 Scale est. = 2.4311 n = 72
What is wrong here? Which one should I proceed with? Is it correct to assume that the differences between the treatments are significant according to the $gam parametric coefficient p-values?
It is a first attempt at additive models trying to catch information from an inter annual fluctuating curve, where a linear model perhaps not completely satisfy.
Observed:
Fitted with the first two models:
The overall goal is still to test for differences between the treatments. Along the entire treatment period, as well as within the years.
For the general model structure, this is the initial approach. Looking at things such as adding by=Treatment or a Year interaction, switching the random Plot factor to same penalty structure, etc.
gam_model3 <- gam(N.P ~ s(as.numeric(Year),bs = "tp", k = 8) + as.numeric(Year)*Treatment +
s(as.numeric(Year), Plot, k=8, bs="fs"), data = leaf, method="REML")
Would be happy for any suggestions from the experienced.
Also looking for the appropriate post-hocs.
EDIT
Reproducible example with modified models:
library(tidyr)
library(dplyr)
library(mgcv)
library(gamm4)
library(performance)
library(ggplot2)
# Generate the data
leaf_rand <- data.frame(
Treatment = rep(factor(c("B", "C", "A", "C", "B", "A", "C", "A", "B")),
times = length(2015:2022)),
Year = rep(2015:2022, each = 9),
Plot = factor(rep(1:9, times = 8)),
N.P = round(c(
15.2, 13.6, 15.8, 14.5, 16.1, 17.5, 19.4, 13.8, 14.5,
16.4, 15.7, 14.1, 14.1, 16.1, 18.4, 19.3, 14.8, 16.9,
17.0, 15.1, 18.9, 17.0, 17.3, 17.0, 19.5, 17.0, 16.0,
18.2, 13.1, 19.7, 17.3, 17.3, 17.7, 22.0, 18.3, 14.9,
17.9, 14.4, 16.3, 15.7, 15.8, 20.0, 14.6, 15.4, 12.8,
17.2, 12.4, 13.9, 11.2, 16.3, 22.0, 18.5, 12.4, 12.1,
16.6, 13.6, 17.0, 10.5, 16.2, 18.3, 17.9, 14.5, 14.0,
22.0, 15.4, 20.0, 17.5, 20.8, 21.8, 21.4, 19.4, 17.1
))
)
### GAMM Models ---------------------------------------------------------------
gam_model <- gamm4(N.P ~ Treatment +
s(Year, k = 5, bs = "tp") +
s(Year, Treatment, k = 5, bs = "fs"),
random = ~(1|Plot),
data = leaf_rand, REML = TRUE)
gam_model1 <- gamm(N.P ~ Treatment +
s(Year, k = 5, bs = "tp") +
s(Year, Treatment, k = 5, bs = "fs"),
random = list(Plot=~1),
data = leaf_rand, REML = TRUE)
gam_model2 <- gam(N.P ~ Treatment +
s(Year, k = 5, bs = "tp") +
s(Year, Treatment, k = 5, bs = "fs") +
s(Plot, k = 9, bs = 're'), data = leaf_rand, method="REML")
# Model summaries
summary(gam_model$gam)
summary(gam_model1$gam)
summary(gam_model2)
# Compare the models
compare <- compare_performance(gam_model$gam, gam_model1$gam, gam_model2)
print(compare)
### ggplot by treatment
ggplot(leaf_rand, aes(x = Year, y = N.P, color = Treatment)) +
geom_point() +
geom_smooth(aes(group = Treatment), method = "loess", se = FALSE) +
theme_minimal() +
theme(legend.position = "bottom") +
labs(
x = "Year",
y = "N.P",
color = "Treatment"
)
What is wrong here?
Which one should I proceed with?
Is it correct to assume that the differences between the treatments are significant according to the $gam parametric coefficient p-values?
It is a first attempt at additive models trying to catch information from an inter annual fluctuating curve, where a linear model perhaps not completely satisfy.
The overall goal is still to test for differences between the treatments. Along the entire treatment period, as well as within the years.
Became Hot Network Question
GAMs and random effects: Significant differences between GAMM and GAMM4 outputs
When introducing a random effect in gamm and gamm4, I receive different p-values for the $gam parametric coefficients. Which seems to be due to a difference in SE calculation(?).
Differences in the smooth term values as well.
Computing a random effect structure in gam gives a comparable result to that of gamm$gam.
The models:
gam_model <- gamm4(N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment,
random = ~(1|Plot), data = leaf, REML = TRUE)
gam_model1 <- gamm(N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment,
random = list(Plot=~1), data = leaf, method = "REML")
gam_model2 <- gam(N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment
+s(Plot, bs="re", k=9), data = leaf, method = "REML")
The outputs:
> summary(gam_model$gam)
Family: gaussian
Link function: identity
Formula:
N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.2879 0.7223 21.167 < 2e-16 ***
TreatmentBCN -0.8325 0.4583 -1.817 0.07389 .
TreatmentControl -1.3039 0.4550 -2.866 0.00561 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(as.numeric(Year)) 4.238 4.238 3.321 0.0133 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.255
lmer.REML = 293.79 Scale est. = 2.4311 n = 72
> summary(gam_model1$gam)
Family: gaussian
Link function: identity
Formula:
N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.2879 1.1649 13.123 <2e-16 ***
TreatmentBCN -0.8325 1.6475 -0.505 0.615
TreatmentControl -1.3039 1.6475 -0.791 0.432
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(as.numeric(Year)) 5.017 5.017 10.39 5.23e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.246
Scale est. = 2.4311 n = 72
> summary(gam_model2)
Family: gaussian
Link function: identity
Formula:
N.P ~ s(as.numeric(Year), bs = "tp", k = 8) + Treatment + s(Plot,
bs = "re", k = 9)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.2879 1.1649 13.123 <2e-16 ***
TreatmentBCN -0.8325 1.6475 -0.505 0.615
TreatmentControl -1.3039 1.6475 -0.791 0.432
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(as.numeric(Year)) 5.017 5.927 8.832 1.3e-06 ***
s(Plot) 5.552 6.000 12.397 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.651 Deviance explained = 71.3%
-REML = 146.9 Scale est. = 2.4311 n = 72
What is wrong here? Which one should I proceed with? Is it correct to assume that the differences between the treatments are significant according to the $gam parametric coefficient p-values?
For more context about the data in question - there are three treatments and three plots per treatment as the repetitions. The treatments have been applied annually during several years and each plot was sampled once per year.
It is a first attempt at additive models trying to catch information from an inter annual fluctuating curve, where a linear model perhaps not completely satisfy.
Observed:
Fitted with the first two models:
The overall goal is still to test for differences between the treatments. Along the entire treatment period, as well as within the years.
For the general model structure, this is the initial approach. Looking at things such as adding by=Treatment or a Year interaction, switching the random Plot factor to same penalty structure, etc.
gam_model3 <- gam(N.P ~ s(as.numeric(Year),bs = "tp", k = 8) + as.numeric(Year)*Treatment +
s(as.numeric(Year), Plot, k=8, bs="fs"), data = leaf, method="REML")
Would be happy for any suggestions from the experienced.
Also looking for the appropriate post-hocs.