Comment puis-je décider de la durée à utiliser dans la régression LOESS dans R?

J'utilise des modèles de régression LOESS en R, et je veux comparer les sorties de 12 modèles différents avec des tailles d'échantillons variables. Je peux décrire les modèles réels plus en détail si cela aide à répondre à la question.

Voici les tailles d'échantillon:

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

Le modèle de régression LOESS est un ajustement de surface, où l'emplacement X et l'emplacement Y de chaque terrain de baseball sont utilisés pour prédire sw, la probabilité de frappe oscillante. Cependant, je voudrais comparer entre les 12 de ces modèles, mais la définition de la même plage (c'est-à-dire la plage = 0,5) produira des résultats différents car il existe une si large gamme de tailles d'échantillon.

Ma question de base est de savoir comment déterminez la portée de votre modèle? Une plage plus élevée lisse davantage l'ajustement, tandis qu'une plage inférieure capture plus de tendances mais introduit un bruit statistique s'il y a trop peu de données. J'utilise une plage plus élevée pour les plus petits échantillons et une plage plus faible pour les plus grands échantillons.

Que devrais-je faire? Quelle est la bonne règle à suivre lors de la définition de l'intervalle pour les modèles de régression LOESS dans R? Merci d'avance!

Une validation croisée est souvent utilisée, par exemple k- fold, si le but est de trouver un ajustement avec le RMSEP le plus bas. Divisez vos données en k groupes et, en laissant chaque groupe à son tour, ajustez un modèle de loess en utilisant les k -1 groupes de données et une valeur choisie du paramètre de lissage, et utilisez ce modèle pour prédire le groupe exclu. Stockez les valeurs prévues pour le groupe exclu, puis répétez jusqu'à ce que chacun des k groupes ait été omis une fois. En utilisant l'ensemble des valeurs prédites, calculez RMSEP. Répétez ensuite le tout pour chaque valeur du paramètre de lissage que vous souhaitez régler. Sélectionnez le paramètre de lissage qui donne le RMSEP le plus bas sous CV.

C'est, comme vous pouvez le voir, assez lourd sur le plan des calculs. Je serais surpris s'il n'y avait pas d'alternative de validation croisée généralisée (GCV) au vrai CV que vous pourriez utiliser avec LOESS - Hastie et al (section 6.2) indiquent que c'est assez simple à faire et couvert dans l'un de leurs exercices .

Je vous suggère de lire les sections 6.1.1, 6.1.2 et 6.2, ainsi que les sections sur la régularisation des splines de lissage (car le contenu s'applique ici aussi) dans le chapitre 5 de Hastie et al. (2009) The Elements of Statistical Learning: Data mining, inference, and prediction . 2e édition. Springer. Le PDF peut être téléchargé gratuitement.

Je suggère de vérifier les modèles additifs généralisés (GAM, voir le package mgcv dans R). J'apprends juste à leur sujet moi-même, mais ils semblent comprendre automatiquement à quel point la «ondulation» est justifiée par les données. Je vois également que vous traitez avec des données binomiales (grève vs pas grève), alors assurez-vous d'analyser les données brutes (c'est-à-dire ne pas agréger en proportions, utiliser les données brutes pas à pas) et utiliser la famille = 'binomial' (en supposant que vous allez utiliser R). Si vous avez des informations sur les pichets et les frappeurs individuels qui contribuent aux données, vous pouvez probablement augmenter votre puissance en faisant un modèle mixte additif généralisé (GAMM, voir le package gamm4 dans R) et en spécifiant le lanceur et le frappeur comme des effets aléatoires (et encore , définissant family = 'binomial'). Finalement, vous voulez probablement permettre une interaction entre les lissées de X & Y, mais je n'ai jamais essayé moi-même, donc je ne sais pas comment faire. Un modèle gamm4 sans l'interaction X * Y ressemblerait à:

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

À bien y penser, vous voudrez probablement laisser les lissages varier à l'intérieur de chaque niveau de type de terrain et de souplesse de la pâte. Cela rend le problème plus difficile car je n'ai pas encore découvert comment laisser les lissages varier selon plusieurs variables de manière à produire ensuite des tests analytiques significatifs ( voir mes requêtes à la liste des modèles mixtes R-SIG ). Tu pourrais essayer:

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)

Mais cela ne donnera pas de tests significatifs des lissages. Pour tenter de résoudre ce problème moi-même, j'ai utilisé un rééchantillonnage bootstrap où, à chaque itération, j'obtiens les prédictions du modèle pour l'espace de données complet, puis calcule les IC bootstap à 95% pour chaque point de l'espace et tous les effets que je souhaite calculer.

Pour une régression de Loess, ma compréhension en tant que non-statisticien, est que vous pouvez choisir votre étendue en fonction de l'interprétation visuelle (un tracé avec de nombreuses valeurs d'étendue peut choisir celui avec le moins de lissage qui semble approprié) ou vous pouvez utiliser la validation croisée (CV) ou validation croisée généralisée (GCV). Ci-dessous se trouve le code que j'ai utilisé pour GCV d'une régression de loess basée sur le code de l'excellent livre de Takezawa, Introduction to Nonparametric Regression (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)
}

et avec mes données, j'ai fait ce qui suit:

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)

Désolé, le code est plutôt bâclé, c'était l'une de mes premières utilisations de R, mais cela devrait vous donner une idée de la façon de faire du GSV pour la régression de Loess pour trouver la meilleure plage à utiliser d'une manière plus objective qu'une simple inspection visuelle. Sur le graphique ci-dessus, vous êtes intéressé par la plage qui minimise la fonction (la plus basse sur la "courbe" tracée).

Si vous passez à un modèle additif généralisé, vous pouvez utiliser la gam()fonction du package mgcv , dans laquelle l'auteur nous assure :

Ainsi, le choix exact de k n'est généralement pas critique: il doit être choisi suffisamment grand pour que vous soyez raisonnablement sûr d'avoir suffisamment de degrés de liberté pour représenter raisonnablement bien la `` vérité '' sous-jacente, mais suffisamment petit pour maintenir une efficacité de calcul raisonnable. Il est clair que «grand» et «petit» dépendent du problème particulier à résoudre.

( kvoici le paramètre des degrés de liberté pour le lissé, qui s'apparente au paramètre de lissage de loess)

Vous pouvez écrire votre propre boucle de validation croisée à partir de zéro en utilisant la loess()fonction du statspackage.

1. Configurez une trame de données de jouets.

   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. Configurez des variables utiles pour gérer la boucle de validation croisée.

   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. Exécutez une forboucle imbriquée itérant sur chaque possibilité de travée span.seqet chaque repli 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. $$CV_{(10)} = \frac{1}{10} \sum_{i=1}^{10} MSE_i$$

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

5. $MSE$

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

6. Tracez vos résultats.

   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)

Utilisez le package locfit. C'est une version légèrement modifiée du loess mais beaucoup plus rapide. Il a également une fonction intégrée pour calculer le gcv http://www.statistik.lmu.de/~leiten/Lehre/Material/GLM_0708/Tutorium/locfit.pdf

Le package fANCOVA fournit un moyen automatisé de calculer la plage idéale en utilisant gcv ou 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))