Интерпретация результатов сплайна

Я пытаюсь подобрать сплайн для GLM с использованием R. После того, как я подгоню сплайн, я хочу иметь возможность взять свою результирующую модель и создать файл моделирования в книге Excel.

Например, допустим, у меня есть набор данных, где y - случайная функция от x, и наклон резко меняется в определенной точке (в данном случае @ 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)

Теперь я подхожу к этому, используя

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

и мои результаты показывают

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

На данный момент, я могу использовать функцию предикторов в r и получить совершенно приемлемые ответы. Проблема в том, что я хочу использовать результаты модели для создания рабочей книги в Excel.

Насколько я понимаю, функция предикта состоит в том, что при заданном новом значении «x» r добавляет этот новый x в соответствующую функцию сплайна (либо в функцию для значений выше 500, либо в функцию для значений ниже 500), а затем принимает этот результат и умножает это с соответствующим коэффициентом и с этой точки зрения обрабатывает это как любой другой модельный термин. Как мне получить эти сплайн-функции?

(Примечание: я понимаю, что гамма-GLM с привязкой к журналу может не подходить для предоставленного набора данных. Я не спрашиваю о том, как и когда подходить для GLM. Я предоставляю этот набор в качестве примера для целей воспроизводимости.)

Вы можете перепроектировать сплайн-формулы без необходимости углубляться в Rкод. Достаточно знать, что

• Сплайн является кусочно-полиномиальной функцией.

• Полиномы степени определяются их значениями в точках .д + $d$$d+1$

• Коэффициенты полинома могут быть получены с помощью линейной регрессии.

Таким образом, вам нужно всего лишь создать точку, разнесенную между каждой парой последовательных узлов (включая неявные конечные точки диапазона данных), предсказать значения сплайнов и регрессировать прогноз по степеням от до . Там будет отдельная формула для каждого базового элемента сплайна в каждом таком узле "корзина". Например, в приведенном ниже примере используются три внутренних узла (для четырех ячеек с узлами) и кубические сплайны ( ), в результате чего получается кубических полиномов, каждый с коэффициентами. Потому что относительно высокие степениx 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$участвуют, обязательно сохранить всю точность в коэффициентах. Как вы можете себе представить, полная формула для любого базового элемента сплайна может быть довольно длинной!

Как я упоминал довольно давно , возможность использовать выходные данные одной программы в качестве входных данных для другой (без ручного вмешательства, которое может привести к невоспроизводимым ошибкам) ​​является полезным навыком статистической коммуникации. Этот вопрос представляет собой хороший пример того, как применяется этот принцип: вместо того, чтобы копировать эти шестнадцатизначных коэффициента вручную, мы можем объединить способ преобразования вычисленных сплайнов в формулы, понятные для Excel. Все, что нам нужно сделать, это извлечь сплайн-коэффициенты, как описано выше, переформатировать их в Excel-подобные формулы, а затем скопировать и вставить их в Excel.$64$RR

Этот метод будет работать с любым статистическим программным обеспечением, даже недокументированным проприетарным программным обеспечением, исходный код которого недоступен.

Вот пример, взятый из вопроса, но модифицированный, чтобы иметь узлы в трех внутренних точках ( ), а также в конечных точках . На графиках показана версия с последующим рендерингом в Excel. Очень мало настроек было выполнено в любой среде (кроме указания цветов, чтобы приблизительно соответствовать цветам Excel по умолчанию).( 1 , 1000 $200, 500, 800$$(1, 1000)$RR

(Вертикальные серые линии сетки в Rверсии показывают, где находятся внутренние узлы.)

Вот полный Rкод. Это несложный хак, полностью полагающийся на pasteфункцию, выполняющую манипуляции со строками. (Лучше было бы создать шаблон формулы и заполнить его с помощью команд сопоставления строк и подстановки.)

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

Первая формула сплайн-вывода (из четырех произведенных здесь)

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

Чтобы это работало в Excel, все, что вам нужно сделать, это удалить окружающие кавычки и поставить перед ними знак «=». (Приложив немного больше усилий, вы могли бы Rнаписать файл, который при импорте в Excel будет содержать копии этих формул во всех нужных местах.) Вставьте его в поле формулы и затем перетаскивайте эту ячейку вокруг, пока «A1» не будет ссылаться на первый значение, где сплайн должен быть вычислен. Скопируйте и вставьте (или перетащите) эту ячейку, чтобы вычислить значения для других ячеек. Я заполнил ячейки B2: E: 102 этими формулами, ссылаясь значения в ячейках A2: A102.$x$$x$

Вы уже сделали следующее:

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

Теперь я покажу вам, как предсказать (ответ) для x = 12 двумя различными способами: во-первых, использовать функцию предиката (самый простой способ!)

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

2-й способ основан на модели матрицы напрямую. Примечание, которое я использовал, expтак как используемая функция ссылки - log.

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

Обратите внимание, что выше я извлек 12-й элемент, поскольку он соответствует x = 12. Если вы хотите сделать прогноз для х вне обучающего набора, то вы просто можете снова использовать функцию прогнозирования. Допустим, мы хотим найти прогнозируемое значение ответа для x = 1100, тогда

> predict(spline1, newdata=data.frame(x=1100),type="response")
       1 
366.3483 

Возможно, вам будет проще использовать усеченную степень мощности для сплайнов кубической регрессии, используя rmsпакет R. Как только вы подгоняете модель, вы можете получить алгебраическое представление подогнанной сплайновой функции, используя функции Functionили .latexrms