Interpretando resultados do spline

20

Estou tentando ajustar um spline para um GLM usando R. Depois de ajustar o spline, quero poder pegar meu modelo resultante e criar um arquivo de modelagem em uma pasta de trabalho do Excel.

Por exemplo, digamos que eu tenho um conjunto de dados em que y é uma função aleatória de xe a inclinação muda abruptamente em um ponto específico (neste caso @ 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)

Agora eu me encaixo nisso usando

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

e meus resultados mostram

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

Neste ponto, eu posso usar a função de previsão dentro de r e obter respostas perfeitamente aceitáveis. O problema é que eu quero usar os resultados do modelo para criar uma pasta de trabalho no Excel.

Meu entendimento da função de previsão é que, dado um novo valor "x", r conecta esse novo x na função spline apropriada (seja a função para valores acima de 500 ou a para valores abaixo de 500), então pega esse resultado e multiplica pelo coeficiente apropriado e a partir desse ponto o trata como qualquer outro termo modelo. Como obtenho essas funções spline?

(Observação: percebo que um GLM gama vinculado a log pode não ser apropriado para o conjunto de dados fornecido. Não estou perguntando como ou quando ajustar os GLMs. Estou fornecendo esse conjunto como um exemplo para fins de reprodutibilidade.)

Eric
fonte
7
Eu sugeriria, se possível, evitar a inclusão de código que exclua todas as variáveis ​​( rm(list=ls())), especialmente não sem nenhum aviso. Alguém pode copiar e colar o código em uma sessão aberta de R onde eles têm algumas variáveis já (mas nenhum chamado x, y, dfou spline1) e falta de que seu código apaga o seu trabalho. É meio idiota para eles fazerem isso? Sim. Mas ainda é educado deixá-los decidir quando excluir suas próprias variáveis.
Glen_b -Reinstate Monica

Respostas:

25

Você pode fazer engenharia reversa das fórmulas de spline sem precisar entrar no Rcódigo. Basta saber que

  • Um spline é uma função polinomial composta por partes.

  • dd+1

  • Os coeficientes de um polinômio podem ser obtidos via regressão linear.

d+1xxdd=34×4=16d+1=4x

64RR

Este método funcionará com qualquer software estatístico, mesmo software proprietário não documentado, cujo código fonte não esteja disponível.

200,500,800(1,1000)RR

Gráficos R

Gráficos em Excel

(As linhas de grade cinzas verticais na Rversão mostram onde estão os nós internos.)


Aqui está o Rcódigo completo . É um hack não sofisticado, confiando inteiramente na pastefunção para realizar a manipulação de strings. (Uma maneira melhor seria criar um modelo de fórmula e preenchê-lo usando comandos de correspondência e substituição de string.)

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

A primeira fórmula de saída de spline (das quatro produzidas aqui) é

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

Rxx

Snippet do Excel

whuber
fonte
2
ns.formula.. você pensa em R ?! Sério, seu método parece muito útil, mas parece irônico ter que invadir um hack para obter esses parâmetros. Seria muito útil para a saída de uma tabela ..
geotheory
Esta pode ser uma pergunta estúpida: mas são 4 splines que você está plotando ou 4 bases de um spline?
Erosennin 28/01
@Erosennin I depende do que você quer dizer com "um spline". Essas quatro curvas são a base de um spline que é cúbico por partes em quatro intervalos e continuamente diferenciável nos três pontos em que esses intervalos se encontram, conforme descrito pelos três pontos de bala que apresentam minha resposta.
whuber
Obrigado! Eu não pretendia dar gorjetas, parece que existem quatro splines (da resposta) e não quatro curvas que são a base. Mais uma vez, estou aqui apenas tentando entender ...
Erosennin 28/01
1
@Erosennin Sem problemas. Talvez isso ajude: o "spline" é qualquer combinação linear dessas quatro curvas que é determinada pelo processo de ajuste de regressão. Outra maneira de dizer: o spline consiste em um espaço vetorial de curvas que pode ser criado usando combinações lineares dessas quatro curvas.
whuber
4

Você já fez o seguinte:

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

Agora vou mostrar como prever (a resposta) para x = 12 de duas maneiras diferentes: Primeiro, usando a função de previsão (a maneira mais fácil!)

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

A segunda maneira é baseada diretamente na matriz do modelo. Nota: usei, exppois a função de link usada é log.

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

Observe que, acima, extraí o 12º elemento, pois corresponde a x = 12. Se você deseja prever um x fora do conjunto de treinamento, basta usar novamente a função de previsão. Digamos que queremos encontrar o valor de resposta previsto para x = 1100 e, em seguida,

> predict(spline1, newdata=data.frame(x=1100),type="response")
       1 
366.3483 
Stat
fonte
Obrigado pela sua resposta! Mas ainda estou confuso: /. Não sei ao certo o que fazer com essa matriz. Por exemplo, se eu tivesse x = 12, a previsão seria y = 68.78721, mas, ao pesquisar 12 nessa matriz, obtive 0,016816392. A interceptação e o coeficiente originais para x <500 são 4,174603 e 3,830416, respectivamente. exp (4,174603 + 3,8304116 * 0,016816392) <> 68,78721. Além disso, como eu obteria valores para x se x não estivesse no conjunto de treinamento?
Eric
Eu mudei minha resposta.
Stat
Eu adicionei um código para o caso em que x não estava no conjunto de treinamento.
Stat
2
Existe uma maneira de obter 366.3483 para x = 1100 sem usar a função de previsão?
Eric
4

Você pode achar mais fácil usar a base de potência truncada para splines de regressão cúbica, usando o rmspacote R. Depois de ajustar o modelo que você pode recuperar a representação algébrica da função spline equipada usando o Functionou latexfunções rms.

Frank Harrell
fonte
Obrigado. Na verdade, li sua resposta aqui stats.stackexchange.com/questions/67607/… antes de postar. Acho que só preciso entender melhor o que posso fazer com o rms.
Eric
A documentação para Function()realmente não diz o que faz. No meu caso (veja detalhes em Rpubs rpubs.com/EmilOWK/rms_splines ), recebo 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>O -2863.7787valor é o primeiro coef no modelo, 245.72672o segundo e o último coef -873.0223não é visto na equação em nenhum lugar. O mesmo se aplica à saída de latex().
Deleet
Functionfunciona Glm()quando você usa rcscomo a função spline. A saída está reformulando o spline da forma mais simples, escrevendo como se as restrições lineares não estivessem lá (mas estão), conforme detalhado nas notas do curso do RMS .
Frank Harrell