Regressão dos mínimos quadrados passo a passo computação em álgebra linear

22

Como prequel para uma pergunta sobre modelos lineares mistos em R, e para compartilhar como referência para os aficionados por estatísticas iniciantes / intermediárias, decidi publicar como um "estilo de perguntas e respostas" independente as etapas envolvidas no cálculo "manual" do coeficientes e valores previstos de uma regressão linear simples.

O exemplo é com o conjunto de dados R embutido mtcars, e seria configurado como milhas por galão consumido por um veículo atuando como variável independente, regredido sobre o peso do carro (variável contínua) e o número de cilindros como um fator com três níveis (4, 6 ou 8) sem interações.

EDIT: Se você está interessado nesta pergunta, você definitivamente encontrará uma resposta detalhada e satisfatória neste post por Matthew Drury fora da CV .

Antoni Parellada
fonte
Quando você diz "computação manual", o que você procura? É relativamente simples mostrar uma série de etapas relativamente simples para obter estimativas de parâmetros e assim por diante (via ortogonalização de Gram-Schmidt, por exemplo, ou pelos operadores SWEEP), mas não é assim que R faz os cálculos internamente; ele (e a maioria dos outros pacotes de estatísticas) usa a decomposição QR (discutida em várias postagens no site - uma pesquisa na decomposição QR mostra várias postagens, algumas das quais você pode obter valor direto)
Glen_b -Reinstate Monica
Sim. Eu acredito que este era endereços muito bem na resposta por MD I provavelmente deve editar o meu post, talvez enfatizando a abordagem geométrica atrás da minha resposta - espaço de coluna, matriz de projeção ...
Antoni Parellada
Sim! @ Matthew Drury Deseja que eu apague essa linha no OP ou atualize o link?
Antoni Parellada
11
Não tenho certeza se você tem esse link, mas isso está intimamente relacionado, e eu realmente amo a resposta de JM. stats.stackexchange.com/questions/1829/…
Haitao Du

Respostas:

51

Nota : Publiquei uma versão expandida desta resposta no meu site .

Você consideraria postar uma resposta semelhante com o mecanismo R real exposto?

Certo! Abaixo a toca do coelho nós vamos.

A primeira camada é lma interface exposta ao programador R. Você pode procurar a fonte disso apenas digitando lmno console R. A maior parte (como a maioria da maioria dos códigos de nível de produção) é ocupada na verificação de entradas, na configuração de atributos de objetos e no lançamento de erros; mas essa linha se destaca

lm.fit(x, y, offset = offset, singular.ok = singular.ok, 
                ...)

lm.fité outra função R, você pode chamá-lo você mesmo. Embora lmtrabalhe convenientemente com fórmulas e quadro de dados, lm.fitdeseja matrizes, então esse é um nível de abstração removido. Verificando a fonte lm.fit, mais trabalho ocupado e a seguinte linha realmente interessante

z <- .Call(C_Cdqrls, x, y, tol, FALSE)

Agora estamos chegando a algum lugar. .Callé a maneira de R chamar o código C. Existe uma função C, C_Cdqrls na fonte R em algum lugar, e precisamos encontrá-la. Aqui está .

Olhando para a função C, novamente, encontramos principalmente verificação de limites, limpeza de erros e trabalho ocupado. Mas essa linha é diferente

F77_CALL(dqrls)(REAL(qr), &n, &p, REAL(y), &ny, &rtol,
        REAL(coefficients), REAL(residuals), REAL(effects),
        &rank, INTEGER(pivot), REAL(qraux), work);

Então agora estamos no nosso terceiro idioma, R chamou C, que está chamando para o fortran. Aqui está o código fortran .

O primeiro comentário diz tudo

c     dqrfit is a subroutine to compute least squares solutions
c     to the system
c
c     (1)               x * b = y

(curiosamente, parece que o nome dessa rotina foi alterado em algum momento, mas alguém se esqueceu de atualizar o comentário). Então, finalmente chegamos ao ponto em que podemos fazer álgebra linear e resolver o sistema de equações. É nesse tipo de coisa que o fortran é realmente bom, o que explica por que passamos por tantas camadas para chegar até aqui.

O comentário também explica o que o código fará

c     on return
c
c        x      contains the output array from dqrdc2.
c               namely the qr decomposition of x stored in
c               compact form.

QR

A primeira coisa que acontece, e de longe a mais importante, é

call dqrdc2(x,n,n,p,tol,k,qraux,jpvt,work)

Isso chama a função fortran dqrdc2em nossa matriz de entrada x. O que é isso?

 c     dqrfit uses the linpack routines dqrdc and dqrsl.

Finalmente chegamos ao linpack . O Linpack é uma biblioteca de álgebra linear fortran que existe desde os anos 70. A álgebra linear mais séria acaba chegando ao linpack. No nosso caso, estamos usando a função dqrdc2

c     dqrdc2 uses householder transformations to compute the qr
c     factorization of an n by p matrix x.

XX=QRQRQR

XtXβ=XtY

muito facilmente. De fato

XtX=RtQtQR=RtR

então todo o sistema se torna

RtRβ=RtQty

RXtX

Rβ=Qty

Rconstant * beta_n = constantβnβQR

Matthew Drury
fonte
4
Este foi o ensaio mais divertido de matemática / codificação que se pode imaginar. Não sei quase nada sobre codificação, mas o seu "passeio" pelas entranhas de uma função aparentemente inócua do R foi verdadeiramente revelador. Excelente escrita! Desde que "gentilmente" fez o truque ... Você poderia gentilmente considerar este como um desafio relacionado? :-)
Antoni Parellada
6
+1 Eu não tinha visto isso antes, bom resumo. Apenas para adicionar um pouco de informação, caso o @Antoni não esteja familiarizado com as transformações de Household; é essencialmente uma transformação linear que permite zerar uma parte da matriz R que você está tentando alcançar sem estragar as partes com as quais você já lidou (contanto que você a passe na ordem certa), tornando-a ideal para transformar matrizes na forma triangular superior (as rotações de Givens fazem um trabalho semelhante e talvez sejam mais fáceis de visualizar, mas são um pouco mais lentas). À medida que você constrói R, você deve, ao mesmo tempo, construir Q
Glen_b
2
Matthew (+1), sugiro que você comece ou termine sua postagem com um link para o seu artigo muito mais detalhado madrury.github.io/jekyll/update/2016/07/07/20/lm-in-R.html .
Ameba diz Reinstate Monica
3
-1 para desbaste e não descendo para o código da máquina.
S. Kolassa - Restabelece Monica
3
(Desculpe,
estou
8

Os cálculos passo a passo reais em R são belamente descritos na resposta de Matthew Drury neste mesmo tópico. Nesta resposta, quero percorrer o processo de provar a si mesmo que os resultados em R com um exemplo simples podem ser alcançados seguindo a álgebra linear de projeções no espaço da coluna e o conceito de erros perpendiculares (produtos pontuais), ilustrados em postagens diferentes , e bem explicado pelo Dr. Strang em Álgebra linear e suas aplicações , e facilmente acessível aqui .

β

mpg=intercept(cyl=4)+β1weight+D1intercept(cyl=6)+D2intercept(cyl=8)[]

D1D2X

attach(mtcars)    
x1 <- wt

    x2 <- cyl; x2[x2==4] <- 1; x2[!x2==1] <-0

    x3 <- cyl; x3[x3==6] <- 1; x3[!x3==1] <-0

    x4 <- cyl; x4[x4==8] <- 1; x4[!x4==1] <-0

    X <- cbind(x1, x2, x3, x4)
    colnames(X) <-c('wt','4cyl', '6cyl', '8cyl')

head(X)
        wt 4cyl 6cyl 8cyl
[1,] 2.620    0    1    0
[2,] 2.875    0    1    0
[3,] 2.320    1    0    0
[4,] 3.215    0    1    0
[5,] 3.440    0    0    1
[6,] 3.460    0    1    0

[]lm

βProjMatrix=(XTX)1XT[ProjMatrix][y]=[RegrCoefs](XTX)1XTy=β

X_tr_X_inv <- solve(t(X) %*% X)    
Proj_M <- X_tr_X_inv %*% t(X)
Proj_M %*% mpg

          [,1]
wt   -3.205613
4cyl 33.990794
6cyl 29.735212
8cyl 27.919934

A idêntico: coef(lm(mpg ~ wt + as.factor(cyl)-1)).

HatMatrix=X(XTX)1XT

HAT <- X %*% X_tr_X_inv %*% t(X)

y^X(XTX)1XTyy_hat <- HAT %*% mpg

cyl <- as.factor(cyl); OLS <- lm(mpg ~ wt + cyl); predict(OLS):

y_hat <- as.numeric(y_hat)
predicted <- as.numeric(predict(OLS))
all.equal(y_hat,predicted)
[1] TRUE
Antoni Parellada
fonte
11
Em geral, na computação numérica, acredito que é melhor resolver a equação linear em vez de calcular a matriz inversa. Então, acho que beta = solve(t(X) %*% X, t(X) %*% y)é, na prática, mais preciso do que solve(t(X) %*% X) %*% t(X) %*% y.
Matthew Drury
R não faz dessa maneira - ele usa uma decomposição QR. Se você vai descrever o algoritmo usado, em um computador duvido que alguém use o que você mostra.
Reponha Monica - G. Simpson
Não após o algoritmo, apenas tentando entender os fundamentos da álgebra linear.
Antoni Parellada
@AntoniParellada Mesmo nesse caso, ainda acho o pensamento em termos de equações lineares mais esclarecedor em muitas situações.
Matthew Drury
11
Dada a relação periférica desse segmento com os objetivos do nosso site, embora tenha valor na ilustração do uso de Rcálculos importantes, gostaria de sugerir que você considere transformá-lo em uma contribuição para o nosso blog.
whuber