Matriz de suavizado de smooth.spline

Question

Estoy tratando de obtener la matriz más suave de mi smooth.spline encajan. Utilizo los datos de densidad mineral ósea de http://statweb.stanford.edu/~tibs/ElemStatLearn/ .

bone <-read.table("bone.data", header=TRUE)
bmd_age <- smooth.spline(bone$age, bone$spnbmd, all.knots=TRUE, cv=TRUE)
bmd_fit <- predict(bmd_age, sort(bone$age))
df <- bmd_age$df

Para obtener una columna de la matriz de suavizado, puedo sustituir el vector de respuesta (bone$spnbmd) por un vector con un solo 1 y el resto lleno de 0's. Es lo que me recomendó el profesor y lo que encontré en internet https://stat.ethz.ch/pipermail/r-help/2006-June/108471.html .

Así que uso

smooth.matrix = function(x){
  n = length(x);
  sm = matrix(0, n, n);  
  for(i in 1:n){
    y = rep(0, n); y[i]=1;
    sm_i = predict(smooth.spline(x, y, df=df),x)$y;
    sm[,i]= sm_i;
  }
  return(sm)
}

sm <- smooth.matrix(bone$age)

Si la matriz de suavizado es correcta, las dos cantidades siguientes deberían ser iguales (ambos valores ajustados del modelo de spline de suavizado).

fromsm <- sm%*%(bone$spnbmd[order(bone$age)])
fromfit <- bmd_fit$y 

Sin embargo, no lo son. Creo que el problema está en la definición de smooth.matrix donde

sm_i = predict(smooth.spline(x, y, df=df),x)$y;

no utiliza el mismo ajuste de suavizado que en bmd_age. He intentado arreglar el grado de libertad, spar, lambda, cv=FALSE, etc. pero no he tenido suerte hasta ahora. ¿Cómo se puede arreglar?

Answer 1

Tras muchas horas de exploración, esto es lo que he encontrado:

Dado que el algoritmo smooth.spline elige spar en lugar de lambda, sólo es posible (más o menos) fijar spar. Sin embargo, lambda es una función de spar y otra matriz variable. Por lo tanto, arreglar spar no arregla necesariamente lambda. No he encontrado una manera fácil de extraer la matriz de suavizado de smooth.spline. Sin embargo, para el propósito de calcular la varianza, el algoritmo proporcionado en https://stat.ethz.ch/pipermail/r-help/2006-June/108471.html (fix spar en lugar de df) es una estimación cercana de la verdadera matriz de suavizado. La varianza calculada a partir de $SyS^T$ , donde $S$ es la matriz de suavizado estimada, se aproxima bastante a la calculada a partir de la matriz de suavizado correcta.

Otro paquete de R "assist" tiene una función "ssr()" que también realiza una regresión spline suave. No es tan potente como smooth.spline. Pero la función incorporada "hat.ssr()" da la verdadera matriz de suavizado del modelo obtenido de "ssr()".

Answer 2

Las respuestas anteriores sólo aproximan la matriz de suavización. Aquí hay una solución que le dará la matriz de suavizado exacta de la función r smooth.spline(). La clave es reconocer que la matriz de suavizado es sólo una función de los valores de $x$ y el parámetro de penalización $\lambda$ , lo que nos permite suavizar un vector $\tilde{y} = (0,0, ..., 0, 1, 0,...,0)^{T}$ y, por tanto, obtener cada columna de la matriz de suavización.

library(splines)
x = seq(0, 100, by=0.1)
y = x*sin(x) + rnorm(length(x), 0, 0.1)

#use cross-validation to choose best smoothing parameter
spar = seq(0.01, 1, by = 0.01)
cv = rep_len(NA, length(spar))
for(i in 1:length(spar)){
    tempfit = smooth.spline(x, y, spar = spar[i], cv=TRUE, all.knots = TRUE)
    cv[i] = tempfit$cv.crit
}

#use the optimal smoothing parameter to produce a final fit
fit = smooth.spline(x, y, spar = spar[which(cv == min(cv))], cv=TRUE, all.knots = TRUE)

#calculate the smoothing matrix
L = matrix(nrow = length(x), ncol = length(x))
for(j in 1:length(x)){
    yi = rep_len(0, length(x))
    yi[j] = 1
    L[,j] = predict(smooth.spline(x, yi, lambda = fit$lambda, cv=TRUE,
                                  all.knots = TRUE), x)$y
}

La matriz $L$ es la matriz de suavizado resultante.

Answer 3

La respuesta aceptada no es correcta aquí - smooth.matrix está funcionando bien.

La única razón fromsm y fromfit no son iguales en el ejemplo anterior se debe a unos paréntesis mal colocados.

Sustituir fromsm <- sm%*%(bone$spnbmd[order(bone$age)]) con fromsm <- (sm%*%bone$spnbmd)[order(bone$age)] y son los mismos.