Jak zdecydować, jakiego zakresu użyć w regresji LOESS w R?

Korzystam z modeli regresji LOESS w R i chcę porównać wyniki 12 różnych modeli o różnych wielkościach próbek. Potrafię opisać rzeczywiste modele bardziej szczegółowo, jeśli pomoże to w udzieleniu odpowiedzi na pytanie.

Oto przykładowe rozmiary:

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

Model regresji LOESS jest dopasowany do powierzchni, gdzie położenie X i położenie Y każdego boiska do baseballu służy do przewidywania sw, prawdopodobieństwa uderzenia wahadłowego. Chciałbym jednak porównać wszystkie 12 modeli, ale ustawienie tego samego zakresu (tj. Zakres = 0,5) przyniesie różne wyniki, ponieważ istnieje tak szeroki zakres wielkości próbek.

Moje podstawowe pytanie brzmi: w jaki sposób określasz rozpiętość swojego modelu? Wyższy zakres bardziej wygładza dopasowanie, podczas gdy niższy zakres przechwytuje więcej trendów, ale wprowadza szum statystyczny, jeśli jest za mało danych. Używam wyższej rozpiętości dla mniejszych próbek i niższej rozpiętości dla większych próbek.

Co powinienem zrobić? Jaka jest dobra zasada przy ustawianiu zakresu dla modeli regresji LOESS w R? Z góry dziękuję!

Często stosuje się walidację krzyżową, na przykład k- krotnie, jeśli celem jest znalezienie dopasowania z najniższym RMSEP. Podziel dane na k grup i pozostawiając po kolei każdą grupę, dopasuj model lessowy, używając k -1 grup danych i wybranej wartości parametru wygładzania, i użyj tego modelu do przewidzenia dla grupy pominiętej. Zapisz prognozowane wartości dla grupy pominiętej, a następnie powtarzaj ją, aż każda z k grup zostanie pominięta jeden raz. Korzystając z zestawu przewidywanych wartości, oblicz RMSEP. Następnie powtórz całość dla każdej wartości parametru wygładzania, który chcesz dostroić. Wybierz ten parametr wygładzania, który daje najniższy RMSEP w CV.

Jest to, jak widać, dość ciężkie obliczeniowo. Byłbym zaskoczony, gdyby nie było uogólnionej alternatywy weryfikacji krzyżowej (GCV) alternatywnej do prawdziwego CV, której można by używać z LOESS - Hastie i wsp. (Sekcja 6.2) wskazują, że jest to dość proste do zrobienia i jest omówione w jednym z ich ćwiczeń .

Proponuję przeczytać sekcje 6.1.1, 6.1.2 i 6.2, a także sekcje dotyczące regulacji regularności wygładzania splajnów (jak również tutaj treść) w rozdziale 5 Hastie i in. (2009) Elementy uczenia statystycznego: eksploracja danych, wnioskowanie i przewidywanie . 2. edycja. Skoczek. Plik PDF można pobrać bezpłatnie.

Sugeruję sprawdzenie uogólnionych modeli addytywnych (GAM, patrz pakiet mgcv w R). Sam się o nich uczę, ale zdają się automatycznie ustalać, jak bardzo „perwersyjność” jest uzasadniona danymi. Widzę również, że masz do czynienia z danymi dwumianowymi (strike vs not strike), więc pamiętaj o analizie surowych danych (tj. Nie agreguj do proporcji, użyj surowych danych pitch-by-pitch) i użyj family = „dwumianowy” (przy założeniu, że będziesz używać R). Jeśli masz informacje na temat tego, co poszczególne miotacze i hittery wnoszą do danych, prawdopodobnie możesz zwiększyć swoją moc, wykonując uogólniony model mieszanego dodatku (GAMM, patrz pakiet gamm4 w R) i określając miotacz i hittera jako efekty losowe (i ponownie , ustawienie rodziny = „dwumianowy”). Wreszcie, prawdopodobnie chcesz pozwolić na interakcję między płynami X i Y, ale sam nigdy tego nie próbowałem, więc nie wiem, jak sobie z tym poradzić. Model gamm4 bez interakcji X * Y wyglądałby następująco:

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

Jeśli się nad tym zastanowić, prawdopodobnie chcesz, aby wygładzenia zmieniały się w zależności od rodzaju wysokości i sposobu trzymania ciasta. To sprawia, że ​​problem jest trudniejszy, ponieważ nie dowiedziałem się jeszcze, w jaki sposób pozwolić wygładzeniom zmieniać się w zależności od wielu zmiennych w sposób, który następnie generuje sensowne testy analityczne ( patrz moje zapytania do listy modeli mieszanych R-SIG ). Możesz spróbować:

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)

Ale to nie da sensownych testów wygładzania. Próbując rozwiązać ten problem sam, użyłem resamplingu bootstrap, gdzie na każdej iteracji uzyskuję prognozy modelu dla pełnej przestrzeni danych, a następnie obliczam 95% CI bootstap dla każdego punktu w przestrzeni i wszelkich efektów, które chcę obliczyć.

W przypadku regresji lessowej, moim zdaniem jako statystykę, możesz wybrać zakres w oparciu o interpretację wizualną (wykres z licznymi wartościami rozpiętości może wybrać ten z najmniejszym wygładzeniem, który wydaje się odpowiedni) lub możesz użyć weryfikacji krzyżowej (CV) lub uogólniona walidacja krzyżowa (GCV). Poniżej znajduje się kod, którego użyłem do GCV regresji lessowej na podstawie kodu z doskonałej książki Takezawy, Wprowadzenie do regresji nieparametrycznej (z 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)
}

i na podstawie moich danych wykonałem następujące czynności:

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)

Niestety, kod jest raczej niechlujny, to był jeden z moich pierwszych przypadków używania R, ale powinien dać ci wyobrażenie o tym, jak zrobić GSV dla regresji lessowej, aby znaleźć najlepszy zakres do użycia w bardziej obiektywny sposób niż zwykła kontrola wizualna. Na powyższym wykresie interesuje Cię rozpiętość, która minimalizuje funkcję (najniższa na narysowanej „krzywej”).

Jeśli przełączysz się na uogólniony model addytywny, możesz użyć gam()funkcji z pakietu mgcv , w którym autor zapewnia nas :

Tak więc dokładny wybór k nie jest na ogół krytyczny: należy wybrać, aby był wystarczająco duży, aby mieć wystarczającą swobodę reprezentowania leżącej u podstaw „prawdy”, ale wystarczająco mały, aby utrzymać rozsądną wydajność obliczeniową. Oczywiście „duży” i „mały” zależą od konkretnego problemu, który ma zostać rozwiązany.

( ktutaj jest parametr stopni swobody dla wygładzacza, który jest podobny do parametru gładkości lessa)

Możesz napisać od początku własną pętlę weryfikacji krzyżowej, która korzysta z loess()funkcji z statspakietu.

1. Skonfiguruj zabawkową ramkę danych.

   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. Ustaw przydatne zmienne do obsługi pętli weryfikacji krzyżowej.

   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. Uruchom zagnieżdżoną forpętlę iterującą po każdej możliwej rozpiętości span.seqi po każdym zagięciu 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. Przedstaw swoje wyniki.

   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)

Użyj pakietu locfit. Jest to nieco zmodyfikowana wersja lessa, ale o wiele szybsza. Posiada również wbudowaną funkcję obliczania gcv http://www.statistik.lmu.de/~leiten/Lehre/Material/GLM_0708/Tutorium/locfit.pdf

FANCOVA pakiet zawiera zautomatyzowany sposób obliczyć zakres idealne pomocą GCV lub 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))