Interpretowanie wyników splajnu

20

Usiłuję dopasować splajn dla GLM za pomocą R. Po dopasowaniu splajnu chcę móc wziąć wynikowy model i utworzyć plik modelowania w skoroszycie programu Excel.

Załóżmy na przykład, że mam zestaw danych, w którym y jest losową funkcją x, a nachylenie zmienia się nagle w określonym punkcie (w tym przypadku @ x = 500).

set.seed(1066)
x<- 1:1000
y<- rep(0,1000)

y[1:500]<- pmax(x[1:500]+(runif(500)-.5)*67*500/pmax(x[1:500],100),0.01)
y[501:1000]<-500+x[501:1000]^1.05*(runif(500)-.5)/7.5

df<-as.data.frame(cbind(x,y))

plot(df)

Teraz dopasowuję to za pomocą

library(splines)
spline1 <- glm(y~ns(x,knots=c(500)),data=df,family=Gamma(link="log"))

i moje wyniki pokazują

summary(spline1)

Call:
glm(formula = y ~ ns(x, knots = c(500)), family = Gamma(link = "log"), 
    data = df)

Deviance Residuals: 
     Min       1Q   Median       3Q      Max  
-4.0849  -0.1124  -0.0111   0.0988   1.1346  

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)             4.17460    0.02994  139.43   <2e-16 ***
ns(x, knots = c(500))1  3.83042    0.06700   57.17   <2e-16 ***
ns(x, knots = c(500))2  0.71388    0.03644   19.59   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.1108924)

    Null deviance: 916.12  on 999  degrees of freedom
Residual deviance: 621.29  on 997  degrees of freedom
AIC: 13423

Number of Fisher Scoring iterations: 9

W tym momencie mogę użyć funkcji przewidywania w obrębie r i uzyskać całkowicie akceptowalne odpowiedzi. Problem polega na tym, że chcę użyć wyników modelu do zbudowania skoroszytu w programie Excel.

Rozumiem, że funkcja przewidywania jest taka, że ​​biorąc pod uwagę nową wartość „x”, r wstawia tę nową x do odpowiedniej funkcji splajnu (funkcja dla wartości powyżej 500 lub funkcja dla wartości poniżej 500), a następnie bierze ten wynik i mnoży według odpowiedniego współczynnika i od tego momentu traktuje go jak każdy inny termin modelowy. Jak uzyskać te funkcje splajnu?

(Uwaga: zdaję sobie sprawę, że logm GLM połączony z logiem może nie być odpowiedni dla dostarczonego zestawu danych. Nie pytam o to, jak i kiedy pasować do GLM. Podam ten zestaw jako przykład do celów odtwarzalności).

Eric
źródło
7
Sugeruję, jeśli to możliwe, unikanie dołączania kodu, który usuwa wszystkie zmienne ( rm(list=ls())), szczególnie nie bez ostrzeżenia. Ktoś może skopiować i wkleić kod do otwartej sesji R, gdzie mają już pewne zmienne (ale żaden zwane x, y, dflub spline1) i miss, że kod wyciera swoją pracę. Czy jest to dla nich trochę głupie? Tak. Ale nadal uprzejmie jest pozwolić im zdecydować, kiedy usunąć własne zmienne.
Glen_b

Odpowiedzi:

25

Możesz poddać inżynierii wstecznej formuły splajnu bez konieczności wchodzenia w Rkod. Wystarczy to wiedzieć

  • Splajn jest fragmentaryczną funkcją wielomianową.

  • Wielomiany stopnia są określone przez ich wartości w punkcie .d + 1dd+1

  • Współczynniki wielomianu można uzyskać za pomocą regresji liniowej.

Dlatego musisz tylko utworzyć punkty między każdą parą kolejnych węzłów (w tym niejawne punkty końcowe zakresu danych), przewidzieć wartości splajnu i zresetować przewidywanie względem potęg od do . W każdym takim „bin” węzła będzie osobna formuła dla każdego elementu bazowego splajnu. Na przykład w poniższym przykładzie występują trzy węzły wewnętrzne (dla czterech przedziałów węzłów) i zastosowano splajny sześcienne ( ), co daje wielomianów sześciennych, każdy o współczynniku . Ponieważ stosunkowo wysokie mocex x d d = 3 4 × 4 = 16 d + 1 = 4 xd+1xxdd=34×4=16d+1=4xsą zaangażowane, konieczne jest zachowanie całej precyzji współczynników. Jak można sobie wyobrazić, pełna formuła dla dowolnego elementu bazowego splajnu może trwać dość długo!

Jak wspomniałem już jakiś czas temu , umiejętność korzystania z danych wyjściowych jednego programu jako danych wejściowych innego programu (bez interwencji ręcznej, która może wprowadzić nieodwracalne błędy) jest przydatną umiejętnością komunikacji statystycznej. To pytanie stanowi dobry przykład zastosowania tej zasady: zamiast ręcznego kopiowania tych szesnastocyfrowych współczynników, możemy zhakować razem sposób konwersji splajnów obliczonych przez na formuły zrozumiałe dla Excela. Wszystko, co musimy zrobić, to wyodrębnić współczynniki splajnu z powyższego opisu, sformatować je do formuł podobnych do Excela, a także skopiować i wkleić je do Excela.64RR

Ta metoda będzie działać z dowolnym oprogramowaniem statystycznym, nawet nieudokumentowanym oprogramowaniem prawnie zastrzeżonym, którego kod źródłowy jest niedostępny.

Oto przykład wzięty z pytania, ale zmodyfikowany tak, aby zawierał węzły w trzech punktach wewnętrznych ( ), a także w punktach końcowych . Wykresy pokazują wersję, a następnie renderowanie w programie Excel. Bardzo niewiele dostosowań przeprowadzono w obu środowiskach (poza określaniem kolorów w przybliżeniu, aby pasowały do ​​domyślnych kolorów Excela).( 1 , 1000 )200,500,800(1,1000)RR

R wykresy

Wykresy Excela

(Pionowe szare linie siatki w Rwersji pokazują, gdzie znajdują się wewnętrzne węzły.)


Oto pełny Rkod. Jest to niewyszukany hack, polegający całkowicie na pastefunkcji umożliwiającej manipulację ciągiem. (Lepszym sposobem byłoby utworzenie szablonu formuły i wypełnienie go za pomocą poleceń dopasowywania i zastępowania ciągów).

#
# Create and display a spline basis.
#
x <- 1:1000
n <- ns(x, knots=c(200, 500, 800))

colors <- c("Orange", "Gray", "tomato2", "deepskyblue3")
plot(range(x), range(n), type="n", main="R Version",
     xlab="x", ylab="Spline value")
for (k in attr(n, "knots")) abline(v=k, col="Gray", lty=2)
for (j in 1:ncol(n)) {
  lines(x, n[,j], col=colors[j], lwd=2)
}
#
# Export this basis in Excel-readable format.
#
ns.formula <- function(n, ref="A1") {
  ref.p <- paste("I(", ref, sep="")
  knots <- sort(c(attr(n, "Boundary.knots"), attr(n, "knots")))
  d <- attr(n, "degree")
  f <- sapply(2:length(knots), function(i) {
    s.pre <- paste("IF(AND(", knots[i-1], "<=", ref, ", ", ref, "<", knots[i], "), ", 
                   sep="")
    x <- seq(knots[i-1], knots[i], length.out=d+1)
    y <- predict(n, x)
    apply(y, 2, function(z) {
      s.f <- paste("z ~ x+", paste("I(x", 2:d, sep="^", collapse=")+"), ")", sep="")
      f <- as.formula(s.f)
      b.hat <- coef(lm(f))
      s <- paste(c(b.hat[1], 
            sapply(1:d, function(j) paste(b.hat[j+1], "*", ref, "^", j, sep=""))), 
            collapse=" + ")
      paste(s.pre, s, ", 0)", sep="")
    })
  })
  apply(f, 1, function(s) paste(s, collapse=" + "))
}
ns.formula(n) # Each line of this output is one basis formula: paste into Excel

Pierwsza formuła wyjściowa splajnu (spośród czterech tu wytworzonych) to

"IF(AND(1<=A1, A1<200), -1.26037447288906e-08 + 3.78112341937071e-08*A1^1 + -3.78112341940948e-08*A1^2 + 1.26037447313669e-08*A1^3, 0) + IF(AND(200<=A1, A1<500), 0.278894459758071 + -0.00418337927419299*A1^1 + 2.08792741929417e-05*A1^2 + -2.22580643138594e-08*A1^3, 0) + IF(AND(500<=A1, A1<800), -5.28222778473101 + 0.0291833541927414*A1^1 + -4.58541927409268e-05*A1^2 + 2.22309136420529e-08*A1^3, 0) + IF(AND(800<=A1, A1<1000), 12.500000000002 + -0.0375000000000067*A1^1 + 3.75000000000076e-05*A1^2 + -1.25000000000028e-08*A1^3, 0)"

Aby działało to w programie Excel, wystarczy usunąć otaczające znaki cudzysłowu i poprzedzić je znakiem „=”. (Przy odrobinie wysiłku możesz Rnapisać plik, który po zaimportowaniu przez program Excel zawiera kopie tych formuł we wszystkich właściwych miejscach.) Wklej go do pola formuły, a następnie przeciągnij tę komórkę, aż „A1” odwołuje się do pierwszej wartość, gdzie ma być obliczony splajn. Skopiuj i wklej (lub przeciągnij i upuść) tę komórkę, aby obliczyć wartości dla innych komórek. Wypełniłem komórki B2: E: 102 tymi wzorami, odnosząc się do wartości w komórkach A2: A102.xx

Fragment programu Excel

Whuber
źródło
2
ns.formula.. czy myślisz w R ?! Poważnie, choć twoja metoda wygląda bardzo przydatna, ale ironiczne wydaje się zhackowanie hacka, aby uzyskać te parametry. Byłoby bardzo przydatne, aby
wygenerować
To może być głupie pytanie: ale czy planujesz 4 splajny, czy 4 podstawy jednego splajnu?
Erosennin
@Erosennin I zależy od tego, co rozumiesz przez „jeden splajn”. Te cztery krzywe są podstawą splajnu, który jest kawałek sześcienny w czterech interwałach i ciągle drugi rozróżnialny w trzech punktach, w których te interwały się spotykają, jak opisano w trzech punktorach, które wprowadzają moją odpowiedź.
whuber
Dzięki! Nie chciałem się dziwić, wygląda na to, że są cztery splajny (z odpowiedzi), a nie cztery krzywe, które są podstawą. Znów jestem tutaj, próbując zrozumieć ...
Erosennin
1
@Erosennin Nie ma problemu. Może to pomoże: „splajn” to dowolna liniowa kombinacja tych czterech krzywych określona przez proces dopasowania regresji. Innymi słowy: splajn składa się z wektorowej przestrzeni krzywych, którą można utworzyć, przyjmując kombinacje liniowe tych czterech krzywych.
whuber
4

Wykonałeś już następujące czynności:

> rm(list=ls())
> set.seed(1066)
> x<- 1:1000
> y<- rep(0,1000)
> y[1:500]<- pmax(x[1:500]+(runif(500)-.5)*67*500/pmax(x[1:500],100),0.01)
> y[501:1000]<-500+x[501:1000]^1.05*(runif(500)-.5)/7.5
> df<-as.data.frame(cbind(x,y))
> library(splines)
> spline1 <- glm(y~ns(x,knots=c(500)),data=df,family=Gamma(link="log"))
> 

Teraz pokażę ci, jak przewidywać (odpowiedź) dla x = 12 na dwa różne sposoby: Najpierw użyj funkcji przewidywania (prosty sposób!)

> new.dat=data.frame(x=12)
> predict(spline1,new.dat,type="response")
       1 
68.78721 

Drugi sposób oparty jest bezpośrednio na macierzy modelu. Uwaga Użyłem, expponieważ użytą funkcją link jest log.

> m=model.matrix( ~ ns(df$x,knots=c(500))) 
> prd=exp(coefficients(spline1) %*% t(m)) 
> prd[12]
[1] 68.78721

Zauważ, że powyżej wyodrębniłem 12. element, ponieważ odpowiada to x = 12. Jeśli chcesz przewidzieć x poza zestawem treningowym, możesz po prostu ponownie użyć funkcji przewidywania. Powiedzmy, że chcemy znaleźć przewidywaną wartość odpowiedzi dla x = 1100

> predict(spline1, newdata=data.frame(x=1100),type="response")
       1 
366.3483 
Stat
źródło
Dziękuję za odpowiedź! Ale nadal jestem zdezorientowany: /. Nie jestem pewien, czy wiem, co zrobić z tą matrycą. Na przykład, jeśli miałem x = 12, to przewidywanie mówi y = 68,78721, ale patrząc w górę 12 z tej macierzy otrzymuję 0,016816392. Pierwotny punkt przecięcia i współczynnik dla x <500 wynoszą odpowiednio 4,174603 i 3,830416. exp (4,174603 + 3,8304116 * 0,016816392) <> 68,78721. Dodatkowo, w jaki sposób uzyskałbym wartości dla x, gdyby x nie było w zestawie treningowym?
Eric
Zmieniłem odpowiedź.
Stat
Dodałem kod dla przypadku, gdy x nie było w zestawie treningowym.
Stat
2
Czy istnieje sposób na uzyskanie 366,3483 dla x = 1100 bez użycia funkcji przewidywania?
Eric
4

Łatwiejsze może być użycie skróconej podstawy mocy dla splajnów regresji sześciennej za pomocą rmspakietu R. Po dopasowaniu modelu można uzyskać reprezentację algebraiczną dopasowanej funkcji splajnu za pomocą funkcji Functionlub latexw rms.

Frank Harrell
źródło
Dziękuję Ci. Właściwie czytam twoją odpowiedź tutaj stats.stackexchange.com/questions/67607/... przed wysłaniem. Chyba potrzebuję lepszego zrozumienia tego, co mogę zrobić z RMS.
Eric
Dokumentacja Function()naprawdę nie mówi, co robi. W moim przypadku (szczegóły na Rpubs rpubs.com/EmilOWK/rms_splines ), dostaję wartość jest pierwszym COEF w modelu, w drugim, a ostatni COEF nie widać nigdzie w równaniu. To samo dotyczy wyników . function(x = NA) {-2863.7787+245.72672* x-0.1391794*pmax(x-10.9,0)^3+0.27835881*pmax(x-50.5,0)^3-0.1391794*pmax(x-90.1,0)^3 } <environment: 0x556156e80db8>-2863.7787245.72672-873.0223latex()
Deleet
Functiondziała, Glm()gdy używasz rcsjako funkcji splajnu. Wyjście przeformułowuje splajn w najprostszej formie, pisząc tak, jakby liniowe ograniczenia ogona nie były (ale są), jak wyszczególniono w notatkach z kursu RMS .
Frank Harrell