Por que estou recebendo previsões diferentes para expansão polinomial manual e usando a função R `poly`?

10

Por que estou recebendo previsões diferentes para expansão polinomial manual e usando a polyfunção R ?

set.seed(0)
x <- rnorm(10)
y <- runif(10)
plot(x,y,ylim=c(-0.5,1.5))
grid()

# xp is a grid variable for ploting
xp <- seq(-3,3,by=0.01)
x_exp <- data.frame(f1=x,f2=x^2)
fit <- lm(y~.-1,data=x_exp)
xp_exp <- data.frame(f1=xp,f2=xp^2)
yp <- predict(fit,xp_exp)
lines(xp,yp)

# using poly function
fit2 <- lm(y~ poly(x,degree=2) -1)
yp <- predict(fit2,data.frame(x=xp))
lines(xp,yp,col=2)

insira a descrição da imagem aqui

Minha tentativa:

  • Parece ser um problema com interceptação, quando eu encaixo o modelo com interceptação, ou seja, não -1no modelo formula, as duas linhas são as mesmas. Mas por que sem a interceptação as duas linhas são diferentes?

  • Outra "correção" é usar rawexpansão polinomial em vez de polinomial ortogonal. Se mudarmos o código para fit2 = lm(y~ poly(x,degree=2, raw=T) -1), formaremos 2 linhas iguais. Mas por que?

Haitao Du
fonte
4
Isso está fora do tópico da sua pergunta, mas você costuma estar muito aberto a comentários. Ao ler seu código, a primeira coisa que noto é que você usa =e <-para atribuição de forma inconsistente. Eu realmente não faria isso, não é exatamente confuso, mas adiciona muito ruído visual ao seu código, sem nenhum benefício. Você deve optar por um ou outro para usar em seu código pessoal e ficar com ele.
Matthew Drury
obrigado por me ajudar na codificação! questão corrigida. @MatthewDrury
Haitao Du
3
Aleatório ponta de acompanhamento para fazer <-menos de um aborrecimento para digitar: alt+-.
JAD 2/17
@JarkoDubbeldam obrigado pela dica de codificação. Eu amo atalhos de teclado
Haitao Du

Respostas:

12

Como você nota corretamente, a diferença original é que, no primeiro caso, você usa os polinômios "brutos", enquanto no segundo caso, usa os polinômios ortogonais. Portanto, se a lmchamada posterior fosse alterada para: fit3<-lm(y~ poly(x,degree=2, raw = TRUE) -1)obteríamos os mesmos resultados entre fite fit3. A razão pela qual obtemos os mesmos resultados neste caso é "trivial"; nos encaixamos exatamente no mesmo modelo em que nos encaixamos fit<-lm(y~.-1,data=x_exp), sem surpresas.

Pode-se facilmente verificar se as matrizes dos dois modelos são iguais all.equal( model.matrix(fit), model.matrix(fit3) , check.attributes= FALSE) # TRUE).


O mais interessante é o motivo pelo qual você obterá os mesmos gráficos ao usar uma interceptação. A primeira coisa a notar é que, ao ajustar um modelo com uma interceptação

  • No caso fit2, simplesmente movemos as previsões do modelo verticalmente; a forma real da curva é a mesma.

  • Por outro lado, incluir uma interceptação no caso de fitresultados não apenas em uma linha diferente em termos de posicionamento vertical, mas com uma forma totalmente diferente em geral.

Podemos ver facilmente isso simplesmente anexando os seguintes ajustes no gráfico existente.

fit_b<-lm(y~. ,data=x_exp)
yp=predict(fit_b,xp_exp)
lines(xp,yp, col='green', lwd = 2)

fit2_b<-lm(y~ poly(x,degree=2, raw = FALSE) )
yp=predict(fit2_b,data.frame(x=xp))
lines(xp,yp,col='blue')

insira a descrição da imagem aqui

OK ... Por que os ajustes de não interceptação são diferentes, enquanto os ajustes de inclusão de interceptação são iguais? A captura está novamente na condição de ortogonalidade.

No caso da fit_bmatriz modelo utilizada conter elementos não ortogonais, a matriz Gram crossprod( model.matrix(fit_b) )está longe de ser diagonal; no caso dos fit2_belementos são ortogonais ( crossprod( model.matrix(fit2_b) )é efetivamente diagonal).

fitfit_b XTXfitfit2fit2_b

A questão interessante é por que os fit_be fit2_bsão os mesmos; afinal todas as matrizes de modelo fit_be fit2_bnão são iguais no valor de face . Aqui, apenas precisamos lembrar disso fit_be fit2_bter as mesmas informações. fit2_bé apenas uma combinação linear do fit_bmodo que seus ajustes resultantes serão os mesmos. As diferenças observadas no coeficiente ajustado refletem a recombinação linear dos valores de fit_b, a fim de obtê-los ortogonais. (veja G. Grothendieck responda aqui também para exemplos diferentes.)

usεr11852
fonte
+2,5 obrigado pela ótima resposta. Para o gráfico final, aprendi com o @kjetilb halvorsen: Uma maneira mais abstrata de descrever isso é que o próprio modelo depende apenas de um certo subespaço linear, ou seja, o espaço da coluna definido pela matriz de design. Mas os parâmetros dependem não apenas desse subespaço, mas da base desse subespaço, dados pelas variáveis ​​específicas utilizadas, ou seja, pelas próprias colunas. As previsões do modelo, por exemplo, dependerão apenas do subespaço linear, não da base escolhida.
Haitao Du
espero que você não se importa, eu reformatado um pouco ..
Haitao Du
@ hxd1011: Sem problemas, obrigado por dedicar um tempo para "pentear" um pouco.
usεr11852