解釋樣條結果

May 15, 2014

我正在嘗試使用 R 為 GLM 擬合樣條曲線。擬合樣條曲線後，我希望能夠採用生成的模型並在 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 中使用 predict 函數並獲得完全可以接受的答案。問題是我想使用模型結果在 Excel 中構建工作簿。

我對預測函數的理解是，給定一個新的“x”值，r 將該新 x 插入適當的樣條函數（值高於 500 的函數或值低於 500 的函數），然後它採用該結果並相乘它通過適當的係數，從那時起，它就像任何其他模型項一樣對待。如何獲得這些樣條函數？

（注意：我意識到對數關聯的伽馬 GLM 可能不適合所提供的數據集。我不是在詢問如何或何時擬合 GLM。我提供該集作為可重複性目的的示例。）

您可以對樣條公式進行逆向工程，而無需進入R代碼。 知道這一點就足夠了

樣條是分段多項式函數。

多項式由它們的值決定點。

多項式的係數可以通過線性回歸獲得。

因此，您只需創建在每對連續節點之間間隔的點（包括數據范圍的隱含端點），預測樣條值，並將預測與取決於. 每個這樣的結“bin”中的每個樣條基元素都有一個單獨的公式。例如，在下面的示例中，有三個內部結（對於四個結箱）和三次樣條（) 被使用，導致三次多項式，每個都有係數。因為功率比較大涉及，必須保持係數中的所有精度。正如您可能想像的那樣，任何樣條基本元素的完整公式都會變得很長！

正如我很久以前提到的，能夠將一個程序的輸出用作另一個程序的輸入（無需人工干預，這會引入不可重現的錯誤）是一種有用的統計溝通技巧。這個問題提供了一個很好的例子來說明該原則是如何應用的：而不是複制那些手動十六位數的係數，我們可以拼湊出一種方法，將計算的樣條曲線轉換R為 Excel 可以理解的公式。我們需要做的就是從上述方法中提取樣條係數R，將它們重新格式化為類似 Excel 的公式，然後將它們複製並粘貼到 Excel 中。

這種方法適用於任何統計軟件，甚至是源代碼不可用的無證專有軟件。

這是取自問題的示例，但修改為在三個內部點處有結（) 以及端點. 繪圖顯示R的版本，然後是 Excel 的渲染。在這兩種環境中都很少進行自定義（除了指定顏色R以大致匹配 Excel 的默認顏色）。

（版本中的垂直灰色網格線R顯示了內部結的位置。）

這是完整的R代碼。這是一個簡單的 hack，完全依賴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 中的值。

引用自：https://stats.stackexchange.com/questions/97845

comments powered by Disqus

解釋樣條結果

相關問答