Dear R users,
I would like to ask you for help with coding up. I would like to recreate a
method from a paper of Malcolm et al. (2014) (MALCOLM, I. A., GIBBINS, C.
N., FRYER, R. J., KEAY, J., TETZLAFF, D. & SOULSBY, C. 2014.
The influence of forestry on acidification and recovery: Insights from
long-term hydrochemical and invertebrate data. Ecological Indicators, 37,
Part B, 317-329.) for multi-site long-term hydrochemistry trends analysis
using GAMMs and mgcv package.
The method description says:
"For a chemical determinant (e.g. water pH) a model was ?tted with: (a)
site included as a factor, (b) a common trend (i.e. a temporal (annual)
smoother common across sites), (c) site-speci?c temporal smoothers measuring
departures from the common trend. Each site-speci?c smoother was constrained
to have the same degrees of freedom assuming common temporal drivers. Three
normally distributed random effects were also included: (i) an
autocorrelated year effect (year treated as a factor) common across sites,
(ii) site-speci?c year effects (year treated as a factor) autocorrelated
within sites but independent between sites (with variance and
autocorrelation assumed common across sites), (iii) an independent residual
term weighted by the square root of the annual sample number to account for
changing sampling effort over time; this weighting, one of several tried,
gave residuals with homogenous variance for all determinants.
[...] the models were ?tted with the variance structure held ?xed at that
estimated by REML and with
the smooths estimated by generalized cross-validation. Smooths with zero edf
(i.e. where the term drops out of the model) where also considered in the
optimisation process".
Let's have some example data:
tmpData = as.data.frame(matrix(data = NA, ncol = 3, nrow = 30))
colnames(tmpData) = c("Determinant", "Year", "Site")
tmpData$Determinant =
c(4.1,4.7,4.5,4.6,4.7,5.6,5.2,
5.6,5.9,5.5,4.9,4.8,4.1,4.2,4.2,4.0,3.9,4.7,4.6,4.5,5.5,5.2,5.1,5.2,5.5,5.7,5.5,5.6,5.8,5.9)
tmpData$Year = c(rep(c(1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001,
2002), 3))
tmpData$Site = c(rep("Site A", 10), rep("Site B", 10), rep("Site C", 10))
So, I would imagine this should start like:
library(mgcv)
model = gamm(Determinant ~ factor(Site) + s(Year, bs = "cr") + ... +
random = list(factor(Year) = ~1, ...),
method = "REML", data = tmpData)
But I know only little about capabilities of mixed effect models and how
possibly code up in this case (c) site-speci?c temporal smoothers measuring
departures from
the common trend, and the random effects in this case.
If anyone could supply me with any tip or a hint, I will be very grateful.
Thank you.
Kasia Sawicka
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
