En guise de préquelle à une question sur les modèles mixtes linéaires dans R, et à partager comme référence pour les aficionados de statistiques débutants / intermédiaires, j'ai décidé de publier en tant que "style Q&A" indépendant les étapes impliquées dans le calcul "manuel" du coefficients et valeurs prédites d'une régression linéaire simple.
L'exemple est avec l'ensemble de données intégré R mtcars
, et serait configuré en miles par gallon consommé par un véhicule agissant comme variable indépendante, régressé sur le poids de la voiture (variable continue), et le nombre de cylindres en tant que facteur à trois niveaux (4, 6 ou 8) sans interactions.
EDIT: Si vous êtes intéressé par cette question, vous trouverez certainement une réponse détaillée et satisfaisante dans ce post de Matthew Drury hors CV .
la source
Réponses:
Remarque : J'ai publié une version étendue de cette réponse sur mon site Web .
Sûr! Dans le terrier du lapin, nous allons.
La première couche est
lm
, l'interface exposée au programmeur R. Vous pouvez regarder la source pour cela en tapant simplementlm
sur la console R. La majorité (comme la majorité de la plupart des codes de niveau de production) est occupée à vérifier les entrées, à définir les attributs des objets et à lancer des erreurs; mais cette ligne ressortlm.fit
est une autre fonction R, vous pouvez l'appeler vous-même. Bien qu'illm
fonctionne de manière pratique avec des formules et un bloc de données, illm.fit
souhaite des matrices, c'est donc un niveau d'abstraction supprimé. Vérification de la sourcelm.fit
, travail plus chargé et ligne vraiment intéressante suivanteNous arrivons maintenant quelque part.
.Call
est la façon dont R appelle le code C. Il y a une fonction C, C_Cdqrls quelque part dans la source R, et nous devons la trouver. Ça y est .En regardant la fonction C, encore une fois, nous trouvons principalement la vérification des limites, le nettoyage des erreurs et le travail occupé. Mais cette ligne est différente
Alors maintenant, nous en sommes à notre troisième langage, R a appelé C qui appelle fortran. Voici le code fortran .
Le premier commentaire dit tout
(Fait intéressant, le nom de cette routine a changé à un moment donné, mais quelqu'un a oublié de mettre à jour le commentaire). Donc, nous sommes enfin au point où nous pouvons faire de l'algèbre linéaire et résoudre le système d'équations. C'est le genre de choses pour lesquelles fortran est vraiment bon, ce qui explique pourquoi nous avons traversé tant de couches pour arriver ici.
Le commentaire explique également ce que le code va faire
Fortran va donc résoudre le système en trouvant la décompositionQR
La première chose qui se produit, et de loin la plus importante, est
Cela appelle la fonction fortran
dqrdc2
sur notre matrice d'entréex
. Qu'est-ce que c'est ça?Nous l'avons donc finalement fait pour linpack . Linpack est une bibliothèque d'algèbre linéaire fortran qui existe depuis les années 70. La plus sérieuse des algèbres linéaires trouve finalement son chemin vers linpack. Dans notre cas, nous utilisons la fonction dqrdc2
C'est là que le travail réel est effectué. Il me faudrait une bonne journée complète pour comprendre ce que fait ce code, il est aussi bas que possible. Mais génériquement, nous avons une matrice et nous voulons la factoriser en un produit X = Q R où Q est une matrice orthogonale et R est une matrice triangulaire supérieure. C'est une chose intelligente à faire, car une fois que vous avez Q et R, vous pouvez résoudre les équations linéaires de régressionX X= Q R Q R Q R
très facilement. Effectivement
de sorte que l'ensemble du système devient
mais est triangulaire supérieur et a le même rang que X t X , donc tant que notre problème est bien posé, il est de rang complet, et nous pouvons tout aussi bien résoudre le système réduitR XtX
Mais voici la chose impressionnante. est triangulaire supérieur, donc la dernière équation linéaire ici est juste , donc la résolution de β n est triviale. Vous pouvez ensuite remonter les lignes, une par une, et les remplacer par les β que vous connaissez déjà, en obtenant à chaque fois une équation linéaire simple à résoudre. Donc, une fois que vous avez Q et R , le tout s'effondre vers ce qu'on appelle la substitution à l'envers , ce qui est facile. Vous pouvez en savoir plus à ce sujet ici , où un petit exemple explicite est entièrement élaboré.R βn β Q R
constant * beta_n = constant
la source
Les calculs réels étape par étape dans R sont magnifiquement décrits dans la réponse de Matthew Drury dans ce même fil. Dans cette réponse, je veux parcourir le processus de prouver à soi-même que les résultats dans R avec un exemple simple peuvent être atteints en suivant l'algèbre linéaire des projections sur l'espace de la colonne et le concept d'erreurs perpendiculaires (produit scalaire), illustré dans différents articles , et joliment expliqué par le Dr Strang dans Algèbre linéaire et ses applications , et facilement accessible ici .
lm
À IDENTIQUE
coef(lm(mpg ~ wt + as.factor(cyl)-1))
.y_hat <- HAT %*% mpg
cyl <- as.factor(cyl); OLS <- lm(mpg ~ wt + cyl); predict(OLS)
:la source
beta = solve(t(X) %*% X, t(X) %*% y)
c'est en pratique plus précis quesolve(t(X) %*% X) %*% t(X) %*% y
.R
pour des calculs importants, je voudrais suggérer que vous envisagiez de le transformer en une contribution à notre blog.