¿Cómo decido qué tramo utilizar en la regresión LOESS en R?

Question

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

Aquí están los tamaños de las muestras:

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 lanzamiento de béisbol se utiliza para predecir la probabilidad de swinging strike. Sin embargo, me gustaría comparar entre los 12 modelos, pero si se establece el mismo intervalo (es decir, intervalo = 0,5) se obtendrán resultados diferentes, ya que hay una amplia gama de tamaños de muestra.

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

¿Qué debo hacer? ¿Cuál es una buena regla general a la hora de establecer el span para los modelos de regresión LOESS en R? Gracias de antemano.

Answer 1

A menudo se utiliza una validación cruzada, por ejemplo k -debido, si el objetivo es encontrar un ajuste con el menor RMSEP. Dividir los datos en k grupos y, dejando de lado cada grupo, ajustar un modelo de loess utilizando el k -1 grupos de datos y un valor elegido del parámetro de suavizado, y utilizar ese modelo para predecir para el grupo que queda fuera. Guarde los valores predichos para el grupo que queda fuera y luego repita hasta que cada uno de los k grupos se ha quedado fuera una vez. Utilizando el conjunto de valores predichos, calcule el RMSEP. A continuación, repita todo el proceso para cada valor del parámetro de suavizado que desee afinar. Seleccione el parámetro de suavizado que dé el menor RMSEP bajo CV.

Esto es, como puedes ver, bastante pesado computacionalmente. Me sorprendería que no existiera una alternativa de validación cruzada generalizada (GCV) al CV verdadero que se pudiera utilizar con LOESS - Hastie et al (sección 6.2) indican que es bastante sencillo de hacer y se trata en uno de sus ejercicios.

Le sugiero que lea la sección 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 se aplica aquí también) en el capítulo 5 de Hastie et al. (2009) Los elementos del aprendizaje estadístico: Minería de datos, inferencia y predicción . 2ª edición. Springer. El PDF se puede descargar gratuitamente.

Answer 2

Sugiero consultar los modelos aditivos generalizados (GAM, véase el paquete mgcv en R). Yo mismo estoy aprendiendo sobre ellos, pero parecen calcular automáticamente cuánta "ondulación" está justificada por los datos. También veo que estás tratando con datos binomiales (strike vs no strike), así que asegúrate de analizar los datos crudos (es decir, no agregues a proporciones, usa los datos crudos de lanzamiento por lanzamiento) y usa family='binomial' (asumiendo que vas a usar R). Si tiene información sobre lo que los lanzadores y bateadores individuales están contribuyendo a los datos, probablemente puede aumentar su poder haciendo un modelo mixto aditivo generalizado (GAMM, ver el paquete gamm4 en R) y especificando el lanzador y el bateador como efectos aleatorios (y de nuevo, estableciendo la familia = 'binomial'). Por último, es probable que quiera 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, probablemente quieras dejar que los alisados 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 todavía no he encontrado la manera de dejar que los suavizados varíen por múltiples variables de una manera que posteriormente produzca pruebas analíticas significativas ( ver mis consultas a la lista R-SIG-Mixed-Models ). Puedes intentarlo:

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 alisados. Para intentar resolver este problema, he utilizado el remuestreo bootstrap, en el que en cada iteración obtengo las predicciones del modelo para todo el espacio de datos y luego calculo los IC del 95% para cada punto del espacio y cualquier efecto que quiera calcular.

Answer 3

Para una regresión de loess, mi comprensión como no estadístico, es que usted puede elegir su span basado en la interpretación visual (parcela con numerosos valores de span puede elegir el que tiene la menor cantidad de suavizado que parece apropiado) o puede utilizar la validación cruzada (CV) o validación cruzada generalizada (GCV). A continuación se muestra el código que he utilizado para GCV de una regresión loess basado en el código del excelente libro de Takezawa, Introducción a la regresión no paramétrica (de la 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)

Lamento que el código sea bastante descuidado, esta fue una de mis primeras veces usando R, pero debería darte una idea de cómo hacer GSV para la regresión de loess para encontrar el mejor tramo a utilizar de una manera más objetiva que la simple inspección visual. En el gráfico anterior, usted está interesado en el tramo que minimiza la función (el más bajo en la "curva" trazada).

Answer 4

Si se cambia a un modelo aditivo generlizado, se podría utilizar el gam() de la función mgcv en el que el autor nos asegura :

Por lo tanto, la elección exacta de k no suele ser crítica: debe elegirse de forma 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. Está claro que "grande" y "pequeño" dependen del problema concreto que se aborde.

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

Answer 5

Puedes escribir tu propio bucle de validación cruzada desde cero que utilice el loess() de la función stats paquete.

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. Establezca variables útiles para manejar el bucle 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. Ejecutar una operación anidada for iterando sobre cada posibilidad de tramo en span.seq y cada pliegue en 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 cuadrático medio medio 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 el tramo que dio lugar a la menor $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)