Extracción de la significación de los términos de suavización del juego

Question

Tengo unos datos del curso del tiempo a los que me gustaría ajustar un gam y tener una fácil interpretación, y con ello me refiero a obtener coefficient estimaciones para cada uno de los spline s.

Estos son mis datos:

df <- data.frame(y = c(0.15,0.17,0.07,0.17,0.01,0.15,0.18,0.04,-0.06,-0.08,0,0.03,-0.27,-0.93,0.04,0.12,0.08,0.15,0.04,0.15,0.03,0.09,0.11,0.13,-0.11,-0.32,-0.7,-0.78,0.07,0.04,0.06,0.12,-0.15,0.05,-0.08,0.14,-0.02,-0.14,-0.24,-0.32,-0.78,-0.81,-0.04,-0.25,-0.09,0.02,-0.13,-0.2,-0.04,0,0.02,-0.05,-0.19,-0.37,-0.57,-0.81),
                 log2.time = rep(c(-1, 0, 1, 1.58,2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58, 6.17, 7.39),4))

Pensé en usar el mgcv R package utilizando este modelo:

fit <- mgcv::gam(y ~ s(log2.time), data = df, method = "REML")

Utilizando summary(fit) me sale:

> summary(fit)

Family: gaussian 
Link function: identity 

Formula:
y ~ s(log2.time)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.10893    0.01681  -6.479 3.65e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
               edf Ref.df     F p-value    
s(log2.time) 4.101  5.036 46.36  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.809   Deviance explained = 82.3%
-REML = -27.673  Scale est. = 0.015829  n = 56

Así que, en general, el log2.time El término de suavización es significativamente diferente de cero. Pero como se ha dicho anteriormente, me gustaría obtener un coefficient para cada uno de los spline s, junto con sus errores estándar y una medida de si son significativamente diferentes de cero.

¿Debería usar cubic splines para este fin, con este modelo:

fit <- mgcv::gam(y ~ s(log2.time,bs='cr'), data = df, method = "REML")

¿Y puede el spline coefficient se obtienen de la siguiente manera:

spline.coefs <- coef(fit)[-1]
spline.knots <- fit$smooth[[1]]$xp

for(s in 2:length(spline.knots)){
  spline.coefs[s-1] <- spline.coefs[s-1]/(spline.knots[s]-spline.knots[s-1])
}

Y entonces el coefficient los errores estándar sean:

spline.coefs.se <- sqrt(diag(vcov(fit, unconditional  = TRUE)))[-1]

Y finalmente, ¿pueden los valores p para el spline coefficient ¿se pueden obtener estimaciones?

Answer 1

La diagonal de la matriz de varianza de los parámetros del modelo contiene las varianzas de todos los coeficientes. Si tomas la raíz cuadrada de los elementos de la diagonal obtendrás los errores estándar del intercepto más los coeficientes del spline:

> fit <- mgcv::gam(y ~ s(log2.time,k=3), data = df, method = "REML")
> vcov(fit)
                 (Intercept) s(log2.time).1 s(log2.time).2
(Intercept)     3.122736e-04  -1.059763e-18  -3.707237e-20
s(log2.time).1 -1.059763e-18   3.968441e-03   3.479426e-08
s(log2.time).2 -3.707237e-20   3.479426e-08   3.122736e-04
> sqrt(diag(vcov(fit)))
   (Intercept) s(log2.time).1 s(log2.time).2 
    0.01767127     0.06299556     0.01767127

Si utiliza vcov(fit, unconditional = TRUE) obtendrá una matriz de covarianza (bayesiana) ajustada por haber estimado el parámetro de suavidad, $\lambda$ . Sin esto, la incertidumbre sobre los coeficientes del spline es demasiado pequeña, ya que supone $\lambda$ es conocido.

Los valores indicados por coef(fit) son los parámetros de las dos funciones de base spline de placa fina utilizadas en el modelo. Sin embargo, no son las pendientes de la spline ajustada; la pendiente o primera derivada de la spline varía a lo largo del rango de $x$ . Puedes estimar la pendiente de la spline en cualquier punto utilizando diferencias finitas si lo deseas, y calcular un intervalo de confianza sobre esa pendiente.

Sólo hay un spline aquí, construido a partir de dos funciones de base (bueno, tres dependiendo de cómo se mire, pero una de ellas es la intercepción). El uso que haces de "spline" y "pendiente" es un poco confuso; si quieres aclararlo, actualizaré esta respuesta en caso de que falten elementos o sean irrelevantes.

Tenga en cuenta que su modelo se ajustó en realidad utilizando k = 3 y luego se eliminó la función de base constante, por lo que el alisador tiene un máximo de 2 grados de libertad. Sin embargo, incluso esto se ha penalizado un poco; si no se desea la penalización (selección de suavidad, o estimación de $\lambda$ ), utilice fx = TRUE .

Actualización:

Si desea un ajuste lineal por partes, puede conseguirlo con la base b-spline bs = 'bs' :

fit <- mgcv::gam(y ~ s(log2.time, bs='bs', k = 3, m = 1), data = df, 
                 method = "REML")

donde m = 1 ahora se refiere a la pedir de la spline, no la pena. En este modelo, el punto de articulación del ajuste está en el centro de los datos. La eliminación de la k = 3 resulta en el ajuste por defecto con k = 10 .

Sin embargo, los coeficientes siguen sin ser las pendientes de las secciones lineales. Deberías ser capaz de extraer la derivada de la spline de la base de la b-spline pero no he encontrado la forma de conseguirlo. mgcv en dar esa información todavía. Sin embargo, puedes conseguirlo con diferencias finitas.

library('mgcv')
devtools::install_github('gavinsimpson/gratia')
library('gratia') # for finite differences-based derivatives
library('ggplot2')
library('cowplot')
theme_set(theme_bw())

## fit piecewise-linear model
fit <- gam(y ~ s(log2.time, bs='bs', m = 1), data = df, 
           method = "REML")
## new locations to evaluate everything at
newdf <- with(df, data.frame(log2.time = seq(min(log2.time), 
                                             max(log2.time),
                                             length = 100)))
## extract model coefs
b <- coef(fit)
## get the Xp matrix, the prediction matrix for the spline
Xp <- predict(fit, newdf, type = "lpmatrix")
## evaluate 1st derivative at new locations
d <- fderiv(fit, newdata = newdf)

## stick everything together
pltdf <- data.frame(log2.time  = newdf$log2.time,
                    spline     = Xp[,-1] %*% b[-1],
                    derivative = d[[1]][['s(log2.time)']][['deriv']])

## individual plots
p1 <- ggplot(pltdf, aes(x = log2.time, y = spline)) +
    geom_line()
p2 <- ggplot(pltdf, aes(x = log2.time, y = derivative)) +
    geom_step()
## combine
plot_grid(p1, p2, align = 'hv', axis = 'lrtb', ncol = 1)