Return to Answer

The summaries indicate a difference between the smooth effects of Year between Treatment levels A and C. At least, the p-values are $< .05$,but you might want to consider correcting for multiple testing. Only the gamm4 model indicates a difference in the parametric effect of Treatment levels A and C. But point estimates are all in the same direction.

I do not see the use of comparing fit between the results of the fitted gam, gamm and gamm4 models. They are essentially equivalent models with identical degrees of freedom (hence the warnings), they just used different estimation and/or optimization approaches.

The summaries indicate a difference between the smooth effects of Year between Treatment levels A and C. At least, the p-values are $< .05$. Only the gamm4 model indicates a difference in the parametric effect of Treatment levels A and C. But point estimates are all in the same direction.

I do not see the use of comparing fit between the fitted gam, gamm and gamm4 models. They are essentially equivalent models with identical degrees of freedom (hence the warnings), they just used different estimation and/or optimization approaches.

As pointed out by Ben Bolker in a comment, PQL versus GHQ for Gaussian responses should not make a difference. However, as shown above, results of mgcv::gam, mgcv::gamm and gamm4::gamm4 differ. These differences can be due to different optimizers used, but such differences would generally be small. The standard error differences for the parametric terms and the difference in F-values for the smooth of Year seem substantial. Posting a reproducible example could be helpful.

With 72 observations, and the presented plots of the observed data, you may want to keep the model simple. Sample size is small, differences between treatments seem slight. A factor smooth might be overdoing it, but whether that is the case is perhaps an empirical question; at least a factor smooth takes up less degrees than a by smooth.

I am not sure why you want to have separate smooths per Plot, as it appears the interest is in the effect of Treatment. Consider keeping only the random intercept term for Plot (it seems to explain most variability, going by the F value of the fitted gam), and modeling only interactions of Year and Treatment, either with a 'linear' interaction as you specified in your model formula, or a by or factor smooth.

As pointed out by Ben Bolker in a comment, PQL versus GHQ for Gaussian responses should not make a difference. However, as shown above, results of mgcv::gam, mgcv::gamm and gamm4::gamm4 differ. Differences can be due to different optimizers used, but such differences would generally be small. The standard error differences for the parametric terms seem substantial, I don't know why, could be due to instability due to relatively small sample size.

If you want to test for differences between the treatments, a by smooth might be more appriate (although it takes up more df than a factor smooth, but I honestly don't know how to use those for testing differences).

Assuming you want to take level A for treatment as the reference category, and check whether each of the other two levels differ, I would take the 'ordered factor' approach. This directly allows you to test whether the parametric and smooth terms differ between the different levels of Treatment:

leaf_rand$oTreatment <- ordered(leaf_rand$Treatment)
gam_model <- gamm4(N.P ~ Treatment + s(Year, k = 5) +
                       s(Year, by = oTreatment, k = 5),
                   random = ~(1|Plot),
                   data = leaf_rand, REML = TRUE)

gam_model1 <- gamm(N.P ~ Treatment + s(Year, k = 5) +
                     s(Year, by = oTreatment, k = 5),
                   random = list(Plot=~1),
                   data = leaf_rand, REML = TRUE)

gam_model2 <- gam(N.P ~ Treatment + s(Year, k = 5) +
                        s(Year, by = oTreatment, k = 5) +
                    s(Plot, bs = 're'), data = leaf_rand, 
                   method="REML")

# Model summaries and plots
summary(gam_model$gam)
summary(gam_model1$gam)
summary(gam_model2)

plot(gam_model$gam)   
plot(gam_model1$gam)   
plot(gam_model2)              

The summaries indicate a difference between the smooth effects of Year between Treatment levels A and C. At least, the p-values are $< .05$,but you might want to consider correcting for multiple testing. Only the gamm4 model indicates a difference in the parametric effect of Treatment levels A and C. But point estimates are all in the same direction.

Inspecting the plots (and edf values) shows that the difference between the two differing levels can be described by a negative linear effect over time.

I do not see the use of comparing fit between the results of the fitted gam, gamm and gamm4 models. They are essentially equivalent models with identical degrees of freedom (hence the warnings), they just used different estimation and/or optimization approaches.

The differences are likely due to the different approaches functions gamm4 and gamm use to approximate the likelihood. nlme (and thereby gam and gamm) uses PQL to approximate the integrands. lme4 (and thereby gamm4) uses Gauss-Hermite quadrature.

(RE)ML estimation of GLMMs requires integrating the random effects out of the model likelihood. There is no closed-form solution or ways to solve this analytically, so numerical methods must be used to approximate the integrals.

From the package documentation of function gamm4::gamm4: "gamm4 is based on gamm from package mgcv, but uses lme4 rather than nlme as the underlying fitting engine via a trick due to Fabian Scheipl. gamm4 is more robust numerically than gamm, and by avoiding PQL gives better performance for binary and low mean count data."

Dimitris Rizopoulos gives a great explanation of PQL and the different ways to numerically approximate the integrals: https://stats.stackexchange.com/a/436711/173546

Edit

As pointed out by Ben Bolker in a comment, PQL versus GHQ for Gaussian responses should not make a difference. However, as shown above, results of mgcv::gam, mgcv::gamm and gamm4::gamm4. These differences can be due to different optimizers used, but such differences would generally be small. Posting a reproducible example could be helpful.

With 72 observations, and the presented plots of the observed data, you may want to keep the model simple. Sample size is small, differences between treatments seem slight. A factor smooth might be overdoing it, but whether that is the case is perhaps an empirical question, but it uses up less degrees than a by smooth.

I am not sure why you want to have separate smooths per Plot, as it appears the interest is in the effect of Treatment. Consider keeping only the random intercept term for Plot (it seems to explain most variability, going by the F value of the fitted gam), and modeling only interactions of Year and Treatment, either with a 'linear' interaction as you specified in your model formula, or a by or factor smooth.

The differences are likely due to the different approaches functions gamm4 and gamm use to approximate the likelihood. nlme (and thereby gam and gamm) uses PQL to approximate the integrands. lme4 (and thereby gamm4) uses Gauss-Hermite quadrature.

(RE)ML estimation of GLMMs requires integrating the random effects out of the model likelihood. There is no closed-form solution or ways to solve this analytically, so numerical methods must be used to approximate the integrals.

From the package documentation of function gamm4::gamm4: "gamm4 is based on gamm from package mgcv, but uses lme4 rather than nlme as the underlying fitting engine via a trick due to Fabian Scheipl. gamm4 is more robust numerically than gamm, and by avoiding PQL gives better performance for binary and low mean count data."

Dimitris Rizopoulos gives a great explanation of PQL and the different ways to numerically approximate the integrals: https://stats.stackexchange.com/a/436711/173546

Edit

As pointed out by Ben Bolker in a comment, PQL versus GHQ for Gaussian responses should not make a difference. However, as shown above, results of mgcv::gam, mgcv::gamm and gamm4::gamm4 differ. These differences can be due to different optimizers used, but such differences would generally be small. The standard error differences for the parametric terms and the difference in F-values for the smooth of Year seem substantial. Posting a reproducible example could be helpful.

With 72 observations, and the presented plots of the observed data, you may want to keep the model simple. Sample size is small, differences between treatments seem slight. A factor smooth might be overdoing it, but whether that is the case is perhaps an empirical question; at least a factor smooth takes up less degrees than a by smooth.

I am not sure why you want to have separate smooths per Plot, as it appears the interest is in the effect of Treatment. Consider keeping only the random intercept term for Plot (it seems to explain most variability, going by the F value of the fitted gam), and modeling only interactions of Year and Treatment, either with a 'linear' interaction as you specified in your model formula, or a by or factor smooth.