¿Cómo decido qué lapso usar en la regresión LOESS en R?

Estoy ejecutando modelos de regresión LOESS en R, y quiero comparar los resultados de 12 modelos diferentes con diferentes tamaños de muestra. Puedo describir los modelos reales con más detalles si me ayuda a responder la pregunta.

Aquí están los tamaños de muestra:

Fastballs vs RHH 2008-09: 2002
Fastballs vs LHH 2008-09: 2209
Fastballs vs RHH 2010: 527 
Fastballs vs LHH 2010: 449

Changeups vs RHH 2008-09: 365
Changeups vs LHH 2008-09: 824
Changeups vs RHH 2010: 201
Changeups vs LHH 2010: 330

Curveballs vs RHH 2008-09: 488
Curveballs vs LHH 2008-09: 483
Curveballs vs RHH 2010: 213
Curveballs vs LHH 2010: 162

El modelo de regresión LOESS es un ajuste de superficie, donde la ubicación X y la ubicación Y de cada campo de béisbol se usan para predecir la probabilidad de golpe oscilante. Sin embargo, me gustaría comparar entre los 12 de estos modelos, pero establecer el mismo intervalo (es decir, intervalo = 0,5) arrojará resultados diferentes ya que existe un rango tan amplio de tamaños de muestra.

Mi pregunta básica es ¿cómo determina la duración de su modelo? Un tramo más alto suaviza más el ajuste, mientras que un tramo más bajo captura más tendencias pero introduce ruido estadístico si hay muy pocos datos. Utilizo un intervalo mayor para tamaños de muestra más pequeños y un intervalo inferior para tamaños de muestra más grandes.

¿Qué tengo que hacer? ¿Cuál es una buena regla general cuando se establece el span para modelos de regresión LOESS en R? ¡Gracias por adelantado!

A menudo se usa una validación cruzada, por ejemplo k- pliegue, si el objetivo es encontrar un ajuste con el RMSEP más bajo. Divida sus datos en k grupos y, dejando a cada grupo a su vez, ajuste un modelo de loess usando los grupos de datos k -1 y un valor elegido del parámetro de suavizado, y use ese modelo para predecir el grupo excluido. Almacene los valores pronosticados para el grupo excluido y luego repita hasta que cada uno de los k grupos se haya omitido una vez. Usando el conjunto de valores predichos, calcule RMSEP. Luego repita todo para cada valor del parámetro de suavizado que desea sintonizar. Seleccione el parámetro de suavizado que proporciona el RMSEP más bajo en CV.

Esto es, como puede ver, bastante computacionalmente pesado. Me sorprendería que no hubiera una alternativa de validación cruzada generalizada (GCV) al CV verdadero que pudiera usar con LOESS. Hastie et al (sección 6.2) indican que esto es bastante simple de hacer y está cubierto en uno de sus ejercicios. .

Le sugiero que lea las secciones 6.1.1, 6.1.2 y 6.2, además de las secciones sobre regularización de splines de suavizado (ya que el contenido también se aplica aquí) en el Capítulo 5 de Hastie et al. (2009) Los elementos del aprendizaje estadístico: minería de datos, inferencia y predicción . 2da edición. Saltador. El PDF se puede descargar de forma gratuita.

Sugiero consultar modelos de aditivos generalizados (GAM, consulte el paquete mgcv en R). Solo estoy aprendiendo acerca de ellos, pero parecen darse cuenta automáticamente de la cantidad de "ondulación" justificada por los datos. También veo que se trata de datos binomiales (huelga versus no huelga), así que asegúrese de analizar los datos sin procesar (es decir, no agregue proporciones, use los datos sin procesar paso por paso) y use family = 'binomial' (suponiendo que vas a usar R). Si tiene información sobre qué lanzadores y bateadores individuales están contribuyendo a los datos, probablemente pueda aumentar su poder haciendo un modelo mixto aditivo generalizado (GAMM, vea el paquete gamm4 en R) y especificando el lanzador y el bateador como efectos aleatorios (y nuevamente , estableciendo family = 'binomial'). Finalmente, probablemente desee permitir una interacción entre los suavizados de X e Y, pero nunca lo he intentado yo mismo, así que no sé cómo hacerlo. Un modelo gamm4 sin la interacción X * Y se vería así:

fit = gamm4(
    formula = strike ~ s(X) + s(Y) + pitch_type*batter_handedness + (1|pitcher) + (1|batter)
    , data = my_data
    , family = 'binomial'
)
summary(fit$gam)

Ahora que lo pienso, es probable que desee dejar que los suavizados varíen dentro de cada nivel de tipo de lanzamiento y mano de bateador. Esto hace que el problema sea más difícil ya que aún no he descubierto cómo permitir que los suavizados varíen según las múltiples variables de una manera que posteriormente produzca pruebas analíticas significativas ( vea mis consultas a la lista de modelos mixtos R-SIG ). Tu podrías intentar:

my_data$dummy = factor(paste(my_data$pitch_type,my_data$batter_handedness))
fit = gamm4(
    formula = strike ~ s(X,by=dummy) + s(Y,by=dummy) + pitch_type*batter_handedness + (1|pitcher) + (1|batter)
    , data = my_data
    , family = 'binomial'
)
summary(fit$gam)

Pero esto no dará pruebas significativas de los suavizados. Al intentar resolver este problema yo mismo, utilicé el remuestreo de bootstrap donde en cada iteración obtengo las predicciones del modelo para el espacio de datos completo y luego calculo los IC del 95% de bootstap para cada punto en el espacio y cualquier efecto que me interese calcular.

Para una regresión loess, mi comprensión como no estadístico, es que puede elegir su rango basado en la interpretación visual (el gráfico con numerosos valores de rango puede elegir el que tenga la menor cantidad de suavizado que parezca apropiado) o puede usar la validación cruzada (CV) o validación cruzada generalizada (GCV). A continuación se muestra el código que utilicé para GCV de una regresión loess basada en el código del excelente libro de Takezawa, Introducción a la regresión no paramétrica (de p219).

locv1 <- function(x1, y1, nd, span, ntrial)
{
locvgcv <- function(sp, x1, y1)
{
    nd <- length(x1)

    assign("data1", data.frame(xx1 = x1, yy1 = y1))
    fit.lo <- loess(yy1 ~ xx1, data = data1, span = sp, family = "gaussian", degree = 2, surface = "direct")
    res <- residuals(fit.lo)

    dhat2 <- function(x1, sp)
    {
        nd2 <- length(x1)
        diag1 <- diag(nd2)
        dhat <- rep(0, length = nd2)

        for(jj in 1:nd2){
            y2 <- diag1[, jj]
            assign("data1", data.frame(xx1 = x1, yy1 = y2))
            fit.lo <- loess(yy1 ~ xx1, data = data1, span = sp, family = "gaussian", degree = 2, surface = "direct")
            ey <- fitted.values(fit.lo)
            dhat[jj] <- ey[jj]
            }
            return(dhat)
        }

        dhat <- dhat2(x1, sp)
        trhat <- sum(dhat)
        sse <- sum(res^2)

        cv <- sum((res/(1 - dhat))^2)/nd
        gcv <- sse/(nd * (1 - (trhat/nd))^2)

        return(gcv)
    }

    gcv <- lapply(as.list(span1), locvgcv, x1 = x1, y1 = y1)
    #cvgcv <- unlist(cvgcv)
    #cv <- cvgcv[attr(cvgcv, "names") == "cv"]
    #gcv <- cvgcv[attr(cvgcv, "names") == "gcv"]

    return(gcv)
}

y con mis datos, hice lo siguiente:

nd <- length(Edge2$Distance)
xx <- Edge2$Distance
yy <- lcap

ntrial <- 50
span1 <- seq(from = 0.5, by = 0.01, length = ntrial)

output.lo <- locv1(xx, yy, nd, span1, ntrial)
#cv <- output.lo
gcv <- output.lo

plot(span1, gcv, type = "n", xlab = "span", ylab = "GCV")
points(span1, gcv, pch = 3)
lines(span1, gcv, lwd = 2)
gpcvmin <- seq(along = gcv)[gcv == min(gcv)]
spangcv <- span1[pgcvmin]
gcvmin <- cv[pgcvmin]
points(spangcv, gcvmin, cex = 1, pch = 15)

Lo sentimos, el código es bastante descuidado, esta fue una de mis primeras veces que usé R, pero debería darle una idea de cómo hacer GSV para la regresión loess para encontrar el mejor lapso para usar de una manera más objetiva que la simple inspección visual. En el gráfico anterior, le interesa el intervalo que minimiza la función (el más bajo en la "curva" trazada).

Si cambia a un modelo de aditivo generalizado, puede usar la gam()función del paquete mgcv , en el que el autor nos asegura :

Por lo tanto, la elección exacta de k generalmente no es crítica: debe elegirse para que sea lo suficientemente grande como para estar razonablemente seguro de tener suficientes grados de libertad para representar la 'verdad' subyacente razonablemente bien, pero lo suficientemente pequeña como para mantener una eficiencia computacional razonable. Claramente 'grande' y 'pequeño' dependen del problema particular que se está abordando.

( kaquí está el parámetro de grados de libertad para el suavizador, que es similar al parámetro de suavidad de loess)

Podría escribir su propio bucle de validación cruzada desde cero que utiliza la loess()función del statspaquete.

1. Configurar un marco de datos de juguete.

   set.seed(4)
   x <- rnorm(n = 500)
   y <- (x)^3 + (x - 3)^2 + (x - 8) - 1 + rnorm(n = 500, sd = 0.5)
   plot(x, y)
   df <- data.frame(x, y)
   

2. Configure variables útiles para manejar el ciclo de validación cruzada.

   span.seq <- seq(from = 0.15, to = 0.95, by = 0.05) #explores range of spans
   k <- 10 #number of folds
   set.seed(1) # replicate results
   folds <- sample(x = 1:k, size = length(x), replace = TRUE)
   cv.error.mtrx <- matrix(rep(x = NA, times = k * length(span.seq)), 
                           nrow = length(span.seq), ncol = k)
   

3. Ejecuta un forbucle anidado que itera sobre cada posibilidad de intervalo span.seqy cada pliegue folds.

   for(i in 1:length(span.seq)) {
     for(j in 1:k) {
       loess.fit <- loess(formula = y ~ x, data = df[folds != j, ], span = span.seq[i])
       preds <- predict(object = loess.fit, newdata = df[folds == j, ])
       cv.error.mtrx[i, j] <- mean((df$y[folds == j] - preds)^2, na.rm = TRUE)
       # some predictions result in `NA` because of the `x` ranges in each fold
    }
   }
   

4. Calcule el error medio de la media de validación cruzada de cada uno de los 10 pliegues:

   $$CV_{(10)} = \frac{1}{10} \sum_{i=1}^{10} MSE_i$$

   cv.errors <- rowMeans(cv.error.mtrx)

5. Encuentre qué período resultó en el más bajo .$MSE$

   best.span.i <- which.min(cv.errors)
   best.span.i
   span.seq[best.span.i]
   

6. Traza tus resultados.

   plot(x = span.seq, y = cv.errors, type = "l", main = "CV Plot")
   points(x = span.seq, y = cv.errors, 
          pch = 20, cex = 0.75, col = "blue")
   points(x = span.seq[best.span.i], y = cv.errors[best.span.i], 
          pch = 20, cex = 1, col = "red")
   
   best.loess.fit <- loess(formula = y ~ x, data = df, 
                           span = span.seq[best.span.i])
   
   x.seq <- seq(from = min(x), to = max(x), length = 100)
   
   plot(x = df$x, y = df$y, main = "Best Span Plot")
   lines(x = x.seq, y = predict(object = best.loess.fit, 
                                newdata = data.frame(x = x.seq)), 
         col = "red", lwd = 2)
   

Use el paquete locfit. Es una versión ligeramente modificada del loess pero mucho más rápido. También tiene una función incorporada para calcular gcv http://www.statistik.lmu.de/~leiten/Lehre/Material/GLM_0708/Tutorium/locfit.pdf

El paquete fANCOVA proporciona una forma automatizada de calcular el lapso ideal usando gcv o aic:

FTSE.lo3 <- loess.as(Index, FTSE_close, degree = 1, criterion = c("aicc", "gcv")[2], user.span = NULL, plot = F)
FTSE.lo.predict3 <- predict(FTSE.lo3, data.frame(Index=Index))