Interpretowanie wyników splajnu

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

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 + $d$$d+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 $d+1$$x$$x^d$$d=3$$4\times 4=16$$d+1=4$$x$są 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.$64$RR

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

(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.$x$$x$

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 

Ł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.