Wie entscheide ich, welche Spanne in der LOESS-Regression in R verwendet werden soll?

Ich verwende LOESS-Regressionsmodelle in R und möchte die Ausgaben von 12 verschiedenen Modellen mit unterschiedlichen Stichprobengrößen vergleichen. Ich kann die tatsächlichen Modelle detaillierter beschreiben, wenn dies bei der Beantwortung der Frage hilfreich ist.

Hier sind die Stichprobengrößen:

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

Das LOESS-Regressionsmodell ist eine Oberflächenanpassung, bei der die X-Position und die Y-Position jedes Baseballfelds verwendet werden, um die Wahrscheinlichkeit eines Swinging-Strikes vorherzusagen. Ich möchte jedoch alle 12 dieser Modelle vergleichen, aber die Einstellung derselben Spanne (dh Spanne = 0,5) führt zu unterschiedlichen Ergebnissen, da es einen so großen Bereich an Stichprobengrößen gibt.

Meine grundlegende Frage ist, wie bestimmen Sie die Spannweite Ihres Modells? Eine höhere Spanne glättet die Anpassung stärker, während eine niedrigere Spanne mehr Trends erfasst, bei zu geringen Daten jedoch statistisches Rauschen erzeugt. Ich verwende eine höhere Spanne für kleinere Stichprobengrößen und eine niedrigere Spanne für größere Stichprobengrößen.

Was sollte ich tun? Was ist eine gute Faustregel beim Festlegen der Spanne für LOESS-Regressionsmodelle in R? Danke im Voraus!

Eine Kreuzvalidierung wird häufig verwendet, zum Beispiel k- fach, wenn das Ziel darin besteht, eine Anpassung mit dem niedrigsten RMSEP zu finden. Teilen Sie Ihre Daten in k Gruppen auf und lassen Sie jede Gruppe der Reihe nach aus. Passen Sie ein Lössmodell unter Verwendung der k -1 Datengruppen und eines ausgewählten Werts des Glättungsparameters an und verwenden Sie dieses Modell zur Vorhersage für die ausgelassene Gruppe. Speichern Sie die vorhergesagten Werte für die ausgelassene Gruppe und wiederholen Sie dies, bis jede der k Gruppen einmal ausgelassen wurde. Berechnen Sie den RMSEP anhand des Satzes der vorhergesagten Werte. Wiederholen Sie das Ganze für jeden Wert des Glättungsparameters, den Sie einstellen möchten. Wählen Sie den Glättungsparameter aus, der unter CV den niedrigsten RMSEP ergibt.

Dies ist, wie Sie sehen können, ziemlich rechenintensiv. Ich wäre überrascht, wenn es keine verallgemeinerte Kreuzvalidierungsalternative (GCV) zu einem echten Lebenslauf gäbe, die Sie mit LOESS verwenden könnten. Hastie et al. (Abschnitt 6.2) geben an, dass dies recht einfach ist und in einer ihrer Übungen behandelt wird .

Ich schlage vor, Sie lesen die Abschnitte 6.1.1, 6.1.2 und 6.2 sowie die Abschnitte zur Regularisierung von Glättungssplines (da der Inhalt auch hier gilt) in Kapitel 5 von Hastie et al. (2009) Die Elemente des statistischen Lernens: Data Mining, Inferenz und Vorhersage . 2. Auflage. Springer. Das PDF kann kostenlos heruntergeladen werden.

Ich schlage vor, generalisierte additive Modelle (GAM, siehe das mgcv-Paket in R) zu prüfen. Ich lerne sie nur selbst kennen, aber sie scheinen automatisch herauszufinden, wie viel "Wackeligkeit" durch die Daten gerechtfertigt ist. Ich sehe auch, dass es sich um Binomialdaten handelt (Strike vs. kein Strike). Analysieren Sie daher unbedingt die Rohdaten (dh aggregieren Sie nicht zu Proportionen, verwenden Sie die Rohdaten Tonhöhe für Tonhöhe) und verwenden Sie family = 'binomial' (unter der Annahme, dass Sie R verwenden werden). Wenn Sie Informationen darüber haben, welche einzelnen Pitcher und Hitter zu den Daten beitragen, können Sie wahrscheinlich Ihre Leistung steigern, indem Sie ein generalisiertes additives Mischmodell (GAMM, siehe gamm4-Paket in R) ausführen und Pitcher und Hitter als zufällige Effekte angeben (und erneut) , setze family = 'binomial'). Endlich, Wahrscheinlich möchten Sie eine Interaktion zwischen den Smooths von X & Y zulassen, aber ich habe das noch nie selbst ausprobiert, daher weiß ich nicht, wie ich das anstellen soll. Ein gamm4-Modell ohne die X * Y-Interaktion würde folgendermaßen aussehen:

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

Wenn Sie sich das überlegen, möchten Sie wahrscheinlich, dass die Glättung je nach Tonhöhenart und Schlagfertigkeit variiert. Dies erschwert das Problem, da ich noch nicht herausgefunden habe, wie man die Glättungen um mehrere Variablen variieren lässt, um anschließend aussagekräftige analytische Tests zu erhalten ( siehe meine Anfragen an die R-SIG-Mixed-Models-Liste ). Du könntest es versuchen:

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)

Dies wird jedoch keine aussagekräftigen Tests der Glättungen ergeben. Bei dem Versuch, dieses Problem selbst zu lösen, habe ich Bootstrap-Resampling verwendet, wobei ich bei jeder Iteration die Modellvorhersagen für den gesamten Datenraum erhalte und dann die 95% -CIs des Bootstaps für jeden Punkt im Raum und alle Effekte berechne, die ich berechnen möchte.

Mein Verständnis als Nicht-Statistiker für eine Löß-Regression ist, dass Sie Ihre Spanne anhand der visuellen Interpretation auswählen können (Diagramm mit zahlreichen Spannenwerten kann diejenige mit der geringsten Glättung auswählen, die angemessen erscheint) oder Sie können die Kreuzvalidierung verwenden (CV) oder generalisierte Kreuzvalidierung (GCV). Unten ist Code, den ich für die GCV einer Löß-Regression verwendet habe, basierend auf Code aus Takezawas ausgezeichnetem Buch " Introduction to Nonparametric Regression" (von S. 219).

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)
}

und mit meinen daten habe ich folgendes gemacht:

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)

Es tut mir leid, dass der Code ziemlich schlampig ist. Dies war eines meiner ersten Male, bei dem ich R verwendet habe, aber er sollte Ihnen eine Vorstellung davon geben, wie Sie GSV für die Löß-Regression durchführen können, um die beste Zeitspanne zu finden, die objektiver als eine einfache visuelle Prüfung ist. Im obigen Diagramm interessiert Sie die Spanne, die die Funktion minimiert (niedrigste in der gezeichneten "Kurve").

Wenn Sie zu einem generierten additiven Modell wechseln, können Sie die gam()Funktion aus dem mgcv- Paket verwenden, in dem der Autor uns versichert :

Daher ist die genaue Wahl von k im Allgemeinen nicht kritisch: Sie sollte so groß gewählt werden, dass Sie einigermaßen sicher sind, genügend Freiheitsgrade zu haben, um die zugrunde liegende „Wahrheit“ angemessen gut darzustellen, aber klein genug, um eine angemessene Recheneffizienz aufrechtzuerhalten. Es ist klar, dass "groß" und "klein" vom jeweiligen Problem abhängen.

( kHier ist der Freiheitsgrad-Parameter für den Glatter, der dem Löß-Parameter für die Glätte entspricht.)

Sie können Ihre eigene Kreuzvalidierungsschleife von Grund auf neu schreiben, die die loess()Funktion aus dem statsPaket verwendet.

1. Richten Sie einen Spielzeugdatenrahmen ein.

   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. Richten Sie nützliche Variablen ein, um Kreuzvalidierungsschleifen zu handhaben.

   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. Führen Sie eine verschachtelte forSchleife aus, die über jede Spannenmöglichkeit in span.seqund jede Falte in iteriert 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. Berechnen Sie den durchschnittlichen quadratischen Kreuzvalidierungsfehler aus jeder der 10 Falten:

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

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

5. Finden Sie heraus, welche Spanne zu der niedrigsten .$MSE$

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

6. Zeichnen Sie Ihre Ergebnisse.

   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)

Verwenden Sie das locfit-Paket. Es ist eine leicht modifizierte Version des Lösses, aber viel schneller. Es hat auch eine eingebaute Funktion zur Berechnung der gcv http://www.statistik.lmu.de/~leiten/Lehre/Material/GLM_0708/Tutorium/locfit.pdf

Das fANCOVA- Paket bietet eine automatisierte Möglichkeit, die ideale Spanne mit gcv oder aic zu berechnen:

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))