2
$\begingroup$

I have annual measurements of these 3 IDs and would like to visualise their overall trend (one curve) using R's gamm function in order to control for autocorrelation within each ID.

> head(dat)
  Time ID  Val
1    1 06 0.66
2    2 06 0.81
3    3 06 0.84
4    4 06 1.09
5    5 06 0.79
6    6 06 1.44

library(mgcv)

mod <- gamm(Val ~ s(Time),
                  random = list(ID = ~1), 
                  correlation = corCAR1(form = ~ Time | ID),
                  data = dat)

> summary(mod$gam)

Family: gaussian 
Link function: identity 

Formula:
Val ~ s(Time)

Parametric coefficients:
            Estimate Std. Error t value      Pr(>|t|)    
(Intercept)   0.9853     0.1559   6.321 0.00000000178 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
        edf Ref.df     F p-value  
s(Time)   1      1 4.162  0.0427 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.151   
  Scale est. = 0.56236   n = 194

mod.fit <- data.frame(ID = dat$ID[complete.cases(dat$Val)],
                            Time = dat$Time[complete.cases(dat$Val)],
                            Val = dat$Val[complete.cases(dat$Val)],
                            fit = predict(mod$gam, type = "response"))

ggplot(mod.fit) + 
  geom_line(aes(x = Time, y = Val, colour = ID)) +
  geom_line(aes(x = Time, y = fit, group = ID)) 

However, I am suprised the fit is so "linear"? Am I doing anything wrong?

enter image description here

Data:

dat <- structure(list(Time = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
                               11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 
                               24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 
                               37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 
                               50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 
                               63L, 64L, 65L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 
                               12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 
                               25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 
                               38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 
                               51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 
                               64L, 65L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 
                               13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 
                               26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 
                               39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 
                               52L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 64L, 
                               65L), ID = c("06", "06", "06", "06", "06", "06", "06", "06", 
                                            "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", 
                                            "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", 
                                            "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", 
                                            "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", 
                                            "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", 
                                            "06", "06", "08", "08", "08", "08", "08", "08", "08", "08", "08", 
                                            "08", "08", "08", "08", "08", "08", "08", "08", "08", "08", "08", 
                                            "08", "08", "08", "08", "08", "08", "08", "08", "08", "08", "08", 
                                            "08", "08", "08", "08", "08", "08", "08", "08", "08", "08", "08", 
                                            "08", "08", "08", "08", "08", "08", "08", "08", "08", "08", "08", 
                                            "08", "08", "08", "08", "08", "08", "08", "08", "08", "08", "08", 
                                            "08", "10", "10", "10", "10", "10", "10", "10", "10", "10", "10", 
                                            "10", "10", "10", "10", "10", "10", "10", "10", "10", "10", "10", 
                                            "10", "10", "10", "10", "10", "10", "10", "10", "10", "10", "10", 
                                            "10", "10", "10", "10", "10", "10", "10", "10", "10", "10", "10", 
                                            "10", "10", "10", "10", "10", "10", "10", "10", "10", "10", "10", 
                                            "10", "10", "10", "10", "10", "10", "10", "10", "10", "10", "10"
                               ), Val = c(0.66, 0.81, 0.84, 1.09, 0.79, 1.44, 0.52, 1.79, 1.24, 
                                          0.95, 0.6, 1.32, 1.48, 1.29, 1.83, 1.95, 1.4, 1.16, 0.81, 0.99, 
                                          0.59, 0, 0.36, 0.64, 0.38, 0.32, 0.34, 0.1, 0.46, 0.36, 0.18, 
                                          0.13, 0.33, 0.37, 0.26, 0.13, 0.36, 0.13, 0.52, 0.51, 0.73, 0.71, 
                                          1.27, 1, 1.5, 0.96, 0.44, 1.21, 0.91, 1.13, 0.63, 0.63, 0.58, 
                                          0.64, 0.75, 0.91, 1.11, 0.91, 0.96, 0.95, 0.61, 0.54, 0.62, 0.42, 
                                          0.54, 1.5, 1.75, 1.54, 1.71, 1.27, 1.65, 0.53, 0.58, 0.84, 1.41, 
                                          1.16, 1.45, 0.96, 0.61, 0.52, 0.69, 0.78, 1.06, 0.78, 0.6, 0.84, 
                                          0.53, 0.87, 0.63, 0.53, 0.53, 0.54, 0.31, 0.37, 0.21, 0.24, 0.3, 
                                          0.22, 0.33, 0.24, 0, 0.29, 0, 0.3, 0.29, 0.47, 0.67, 1.13, 1.52, 
                                          2.73, 1.39, 0.53, 1.1, 0.91, 0.94, 0.89, 0.85, 1.33, 0.94, 1.3, 
                                          1.24, 1.33, 1.69, 1.49, 0.76, 0.82, 0.44, 0.63, 0.44, 0.45, 0.74, 
                                          0.72, 1.32, 2.83, 2.16, 3.75, 3.08, 4.2, 4.4, 3.9, 2.35, 3.57, 
                                          2.94, 2.96, 3.28, 2.91, 3.59, 3.38, 2.1, 2.1, 1.4, 1.14, 2.04, 
                                          2.01, 1.13, 1.45, 1.31, 0.92, 0.97, 0.63, 1, 0.91, 0.54, 0.83, 
                                          0.68, 0.96, 0.77, 0.56, 0.74, 1, 1.27, 0.82, 1.85, 1.91, 1.56, 
                                          0.52, 0.27, 0.85, 0.78, 0.85, 0.45, 0.9, 0, 0.64, 1.1, 0.6, 0.42, 
                                          0.25, 0.28, 0.32, 0.57, 0.23, 0.18, 0.14, NA)), class = "data.frame", row.names = c(NA, 
                                                                                                                              -195L))
$\endgroup$
  • 1
    $\begingroup$ quick thoughts: (1) try increasing the maximum dimension (e.g. s(Time,k=30)); (2) the random effect only includes an among-ID shift in the intercept, doesn't allow for different curves by ID; (3) the random effects term is not going to work very well with only three groups - try comparing a simple gam() [with a time-by-ID interaction] $\endgroup$ – Ben Bolker Oct 1 '18 at 12:59
  • 1
    $\begingroup$ I don't think it makes sense to combine a temporal correlation model and a GAM with a time smoother. Currently the wiggliness is modelled by the correlation structure. $\endgroup$ – Roland Oct 1 '18 at 13:17
  • 1
    $\begingroup$ With just three subjects, you could do dat$ID <- factor(dat$ID); mod <- gam(Val ~ ID + s(Time, by = ID), data = dat, select = TRUE). If you had more subjects, a random effects smoother might be sensible: mod <- gam(Val ~ s(Time) + s(ID, Time, bs = "re") + s(ID, bs = "re"), data = dat, select = TRUE). $\endgroup$ – Roland Oct 1 '18 at 13:24
  • 1
    $\begingroup$ @Roland Re your random effects smoother, those s(ID, Time, bs = 're') are really just linear random effects (I guess slopes in this instance). Did you mean to suggest bs = 'fs', which would give random smooths. $\endgroup$ – Gavin Simpson Oct 2 '18 at 2:56
  • 1
    $\begingroup$ Further to @Roland's point about including both a smooth of time and a correlation structure, the issue is that these two components are often unidentifiable; strong autocorrelation about a simple linear trend or a wiggly trend and no autocorrelation are two ways you could describe the data. Unless you have a priori reason to constrain one of the two processes (say by saying the smooth is not very wiggly, or that the correlation parameter is small) then there is often little the model can do to separate these two processes. $\endgroup$ – Gavin Simpson Oct 2 '18 at 2:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy