Sorry about the errors (typos, not syntax errors) - I was forgetting that you'd need to use `gamm()` and hence access the `$gam` component
I don't follow the point about a factor trending up or down. You shouldn't try to use the `$lme` part of the model for this. `summary(mod$gam)` should be sufficient, but as it relates to a spline, this is more a test of whether the spline is different from a horizontal, flat, null line. The problem with splines is that the trend need not just be "trending up or down". In the past, to convey where change in the trend occurs I have used the first derivative of the fitted spline and looked for where in time the 95% confidence interval on the first derivative of the spline doesn't include zero; that shows the regions in time where the trend is significantly increasing or decreasing. I cover how to do this in a blog post I wrote: http://www.fromthebottomoftheheap.net/2011/06/12/additive-modelling-and-the-hadcrut3v-global-mean-temperature-series/ the post contains links to the R code used for the derivatives etc, though it is a little more complex in the case of a model with a trend spline and seasonal spline. I'm supposed to have updated those codes and the post because several people have asked me how I do the analysis for models with multiple spline terms. If you can't get the code to work for your models, ping me back and I'll try to move that to the top of my TO DO list. Note the `Xs(time)Fx1` entries in the `summary(mod$lme)` table refer to the basis functions that represent the spline or at least to some part of those basis functions. You can't really make much practical use out of those values are they relate specifically to way the penalised regression spline model has been converted into an equivalent linear mixed effect form. HTH G On 26 March 2014 12:10, Jacob Cram <cramj...@gmail.com> wrote: > Thanks again Gavin, this works. > gamm() also models the long term trend with a spline s(Time), which is > great. I would still like though, to be able to say whether the factor is > trending up or down over time. Would it be fair to query > summary(mod$lme)$tTable > and to look at the p-value and "Value" corresponding Xs(time)Fx1 value to > identify such a trend? > > Also, here are a few syntax corrections on the code provided in the last > email: > 1)visual appoach > plot(mod$gam, pages = 1) > > 2) quantitative approach > pred <- predict(mod$gam, newdata = newdat, type = "terms") > take <- which(pred[,"s(DoY)"] == max(pred[,"s(DoY)"])) > or > take <- as.numeric(which.max(pred[,"s(DoY)"])) > > Cheers, > -Jacob > > > On Wed, Mar 26, 2014 at 9:46 AM, Gavin Simpson <ucfa...@gmail.com> wrote: >> >> 1) Visually - unless it actually matters exactly on which day in the >> year the peak is observed? If visually is OK, just do `plot(mod, pages >> = 1)` to see the fitted splines on a single page. See `?plot.gam` for >> more details on the plot method. >> >> 2) You could generate some new data to predict upon as follows: >> >> newdat <- data.frame(DoY = seq_len(366), time = mean(foo$time)) >> >> Then predict for those new data but collect the individual >> contributions of the spline terms to the predicted value rather than >> just the final prediction >> >> pred <- predict(mod, newdata = newdat, type = "terms") >> >> Then find the maximal value of the DoY contribution >> >> take <- which(pred$DoY == max(pred$DoY)) >> newdat[take, , drop = FALSE] >> >> You could use >> >> take <- which.max(pred$DoY) >> >> instead of the `which()` I used, but only if there is a single maximal >> value. >> >> This works because the spline terms in the additive model are just >> that; additive. Hence because you haven't let the DoY and time splines >> interact (in the simple model I mentioned, it is more complex if you >> allow these to interact as you then need to predict DoY for each years >> worth of time points), you can separate DoY from the other terms. >> >> None of the above code has been tested, but was written off top of my >> head, but should work or at least get you pretty close to something >> that works. >> >> HTH >> >> G >> >> On 26 March 2014 10:02, Jacob Cram <cramj...@gmail.com> wrote: >> > Thanks Gavin, >> > This seems like a promising approach and a first pass suggests it >> > works >> > with this data. I can't quite figure out how I would go about >> > interrogating >> > the fitted spline to deterine when the peak value happens with respect >> > to >> > DoY. Any suggestions? >> > -Jacob >> > >> > >> > On Tue, Mar 25, 2014 at 9:06 PM, Gavin Simpson <ucfa...@gmail.com> >> > wrote: >> >> >> >> I would probably attack this using a GAM modified to model the >> >> residuals as a stochastic time series process. >> >> >> >> For example >> >> >> >> require("mgcv") >> >> mod <- gamm(y ~ s(DoY, bs = "cc") + s(time), data = foo, >> >> correlation = corCAR1(form = ~ time)) >> >> >> >> where `foo` is your data frame, `DoY` is a variable in the data frame >> >> computed as `as.numeric(strftime(RDate, format = "%j"))` and `time` is >> >> a variable for the passage of time - you could do `as.numeric(RDate)` >> >> but the number of days is probably large as we might encounter more >> >> problems fitting the model. Instead you might do `as.numeric(RDate) / >> >> 1000` say to produce values on a more manageable scale. The `bs = >> >> "cc"` bit specifies a cyclic spline applicable to data measured >> >> throughout a year. You may want to fix the start and end knots to be >> >> days 1 and days 366 respectively, say via `knots = list(DoY = >> >> c(0,366))` as an argument to `gam()` [I think I have this right, >> >> specifying the boundary knots, but let me know if you get an error >> >> about the number of knots]. The residuals are said to follow a >> >> continuois time AR(1), the irregular-spaced counter part to the AR(1), >> >> plus random noise. >> >> >> >> There may be identifiability issues as the `s(time)` and `corCAR1()` >> >> compete to explain the fine-scale variation. If you hit such a case, >> >> you can make an educated guess as to the wiggliness (degrees of >> >> freedom) for the smooth terms based on a plot of the data and fix the >> >> splines at those values via argument `k = x` and `fx = TRUE`, where >> >> `x` in `k = x` is some integer value. Both these go in as arguments to >> >> the `s()` functions. If the trend is not very non-linear you can use a >> >> low value 1-3 here for x and for the DoY term say 3-4 might be >> >> applicable. >> >> >> >> There are other ways to approach this problem of identifiability, but >> >> that would require more time/space here, which I can go into via a >> >> follow-up if needed. >> >> >> >> You can interrogate the fitted splines to see when the peak value of >> >> the `DoY` term is in the year. >> >> >> >> You can also allow the seasonal signal to vary in time with the trend >> >> by allowing the splines to "interact" in a 2d-tensor product spline. >> >> Using `te(DoY, time, bs = c("cc","cr"))` instead of the two `s()` >> >> terms (or using `ti()` terms for the two "marginal" splines and the >> >> 2-d spline). Again you can add in the `k` = c(x,y), fx = TRUE)` to the >> >> `te()` term where `x` and `y` are the dfs for each dimension in the >> >> `te()` term. It is a bit more complex to do this for `ti()` terms. >> >> >> >> Part of the reason to prefer a spline for DoY for the seasonal term is >> >> that one might not expect the seasonal cycle to be a symmetric cycle >> >> as a cos/sin terms would imply. >> >> >> >> A recent ecological paper describing a similar approach (though using >> >> different package in R) is that of Claire Ferguson and colleagues in J >> >> Applied Ecology (2008) http://doi.org/10.1111/j.1365-2664.2007.01428.x >> >> (freely available). >> >> >> >> HTH >> >> >> >> G >> >> >> >> On 25 March 2014 19:14, Jacob Cram <cramj...@gmail.com> wrote: >> >> > Hello all, >> >> > I am thinking about applying season::cosinor() analysis to some >> >> > irregularely spaced time series data. The data are unevenly spaced, >> >> > so >> >> > usual time series methods, as well as the nscosinor() function are >> >> > out. >> >> > My >> >> > data do however trend over time and I am wondering if I can feed date >> >> > as >> >> > a >> >> > variable into my cosinor analyis. In the example below, then I'd >> >> > conclude >> >> > then that the abundances are seasonal, with maximal abundance in mid >> >> > June >> >> > and furthermore, they are generally decreasing over time. >> >> > Can I use both time variables together like this? If not, is there >> >> > some >> >> > better approach I should take? >> >> > Thanks in advance, >> >> > -Jacob >> >> > >> >> > For context lAbundance is logg-odds transformed abundance data of a >> >> > microbial species in a given location over time. RDate is the date >> >> > the >> >> > sample was collected in the r date format. >> >> > >> >> > >> >> >> res <- cosinor(lAbundance ~ RDate, date = "RDate", data = lldata) >> >> > >> >> >>summary(res) >> >> > Cosinor test >> >> > Number of observations = 62 >> >> > Amplitude = 0.58 >> >> > Phase: Month = June , day = 14 >> >> > Low point: Month = December , day = 14 >> >> > Significant seasonality based on adjusted significance level of 0.025 >> >> > = >> >> > TRUE >> >> > >> >> >>summary(res$glm) >> >> > >> >> > >> >> > Call: >> >> > glm(formula = f, family = family, data = data, offset = offset) >> >> > >> >> > Deviance Residuals: >> >> > Min 1Q Median 3Q Max >> >> > -2.3476 -0.6463 0.1519 0.6574 1.9618 >> >> > >> >> > Coefficients: >> >> > Estimate Std. Error t value Pr(>|t|) >> >> > (Intercept) 0.0115693 1.8963070 0.006 0.99515 >> >> > RDate -0.0003203 0.0001393 -2.299 0.02514 * >> >> > cosw -0.5516458 0.1837344 -3.002 0.00395 ** >> >> > sinw 0.1762904 0.1700670 1.037 0.30423 >> >> > --- >> >> > Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 >> >> > >> >> > (Dispersion parameter for gaussian family taken to be 0.9458339) >> >> > >> >> > Null deviance: 70.759 on 61 degrees of freedom >> >> > Residual deviance: 54.858 on 58 degrees of freedom >> >> > AIC: 178.36 >> >> > >> >> > Number of Fisher Scoring iterations: 2 >> >> > >> >> > >> >> > llldata as a csv below >> >> > "RDate","lAbundance" >> >> > 2003-03-12,-3.3330699059335 >> >> > 2003-05-21,-3.04104625745886 >> >> > 2003-06-17,-3.04734680029566 >> >> > 2003-07-02,-4.18791034708572 >> >> > 2003-09-18,-3.04419201802053 >> >> > 2003-10-22,-3.13805060873929 >> >> > 2004-02-19,-3.80688269144794 >> >> > 2004-03-17,-4.50755507726145 >> >> > 2004-04-22,-4.38846502542992 >> >> > 2004-05-19,-3.06618649442674 >> >> > 2004-06-17,-5.20518774876304 >> >> > 2004-07-14,-3.75041853151097 >> >> > 2004-08-25,-3.67882486716196 >> >> > 2004-09-22,-5.22205827512234 >> >> > 2004-10-14,-3.99297508670535 >> >> > 2004-11-17,-4.68793287601157 >> >> > 2004-12-15,-4.31712380781011 >> >> > 2005-02-16,-4.30893550479904 >> >> > 2005-03-16,-4.05781773988454 >> >> > 2005-05-11,-3.94746237402035 >> >> > 2005-07-19,-4.91195185391358 >> >> > 2005-08-17,-4.93590576323119 >> >> > 2005-09-15,-4.85820800095518 >> >> > 2005-10-20,-5.22956391101343 >> >> > 2005-12-13,-5.12244047315448 >> >> > 2006-01-18,-3.04854660925046 >> >> > 2006-02-22,-6.77145858348375 >> >> > 2006-03-29,-4.33151493849021 >> >> > 2006-04-19,-3.36152357710535 >> >> > 2006-06-20,-3.09071584142593 >> >> > 2006-07-25,-3.31430484483825 >> >> > 2006-08-24,-3.09974933041469 >> >> > 2006-09-13,-3.33288992218458 >> >> > 2007-12-17,-4.19942661980677 >> >> > 2008-03-19,-3.86146499633625 >> >> > 2008-04-22,-3.36161599919095 >> >> > 2008-05-14,-4.30878307213324 >> >> > 2008-06-18,-3.74372448768828 >> >> > 2008-07-09,-4.65951429661651 >> >> > 2008-08-20,-5.35984647704619 >> >> > 2008-09-22,-4.78481898261137 >> >> > 2008-10-20,-3.58588161980965 >> >> > 2008-11-20,-3.10625125552057 >> >> > 2009-02-18,-6.90675477864855 >> >> > 2009-03-11,-3.43446932013368 >> >> > 2009-04-23,-3.82688066341466 >> >> > 2009-05-13,-4.44885332005661 >> >> > 2009-06-18,-3.97671552612412 >> >> > 2009-07-09,-3.40185954003936 >> >> > 2009-08-19,-3.44958231694091 >> >> > 2009-09-24,-3.86508161094726 >> >> > 2010-01-28,-4.95281587967569 >> >> > 2010-02-11,-3.78064756876257 >> >> > 2010-03-24,-3.5823501176064 >> >> > 2010-04-27,-4.33635715879999 >> >> > 2010-05-17,-3.90545735473055 >> >> > 2010-07-21,-3.3147176517321 >> >> > 2010-08-11,-4.53218360860017 >> >> > 2010-10-21,-6.90675477864855 >> >> > 2010-11-23,-6.90675477864855 >> >> > 2010-12-16,-6.75158176094352 >> >> > 2011-01-11,-6.90675477864855 >> >> > >> >> > [[alternative HTML version deleted]] >> >> > >> >> > _______________________________________________ >> >> > R-sig-ecology mailing list >> >> > R-sig-ecology@r-project.org >> >> > https://stat.ethz.ch/mailman/listinfo/r-sig-ecology >> >> >> >> >> >> >> >> -- >> >> Gavin Simpson, PhD >> > >> > >> >> >> >> -- >> Gavin Simpson, PhD > > -- Gavin Simpson, PhD _______________________________________________ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
