Calcul matricial pentru estimarea parametrilor modelului de regresie liniară


Dacă, în regresia simplă, pentru estimarea parametrilor β̂₀ și β̂₁, și a erorilor lor standard, calculele sunt simple și ușor de memorat, situația este mult mai complexă la trecerea la regresie liniară multiplă, și devine rapid complicată, pe măsură ce se adaugă predictori. Totuși, se pot simplifica mult și memora procedee de estimare, dacă se folosește calculul matricial. Din fericire pentru studentul testat cu pixul pe foaie, foarte puține astfel de calcule se pot face în timpul scurt dedicat unui examen.

Estimarea parametrilor modelului de regresie multiplă folosind matrici

Folosind exemplul din celelalte capitole, să spunem că dorim să modelăm, pentru a previziona, PIBul, folosind ca predictori Populația, Procentul participării la forța de muncă și Procentul de absolvenți de studii universitare. Acest capitol folosește un eșantion selectat aleator la fiecare 30 de minute.

Modelul

  • \( PIB_{it} = \beta_0 + \beta_1 \times Pop_{it} + \beta_2 \times PPFM_{it} + \beta_3 \times PASU_{it} + \varepsilon_{it} \)
  • \( PIB_{ro} = \beta_0 + \beta_1 \times Pop_{ro} + \beta_2 \times PPFM_{ro} + \beta_3 \times PASU_{ro} + \varepsilon_{ro} \)
  • \( PIB_{es} = \beta_0 + \beta_1 \times Pop_{es} + \beta_2 \times PPFM_{es} + \beta_3 \times PASU_{es} + \varepsilon_{es} \)
  • etc.

La modul general, notând variabila previzionată cu Y și predictorii cu X₁, X₂ etc, vom scrie

  • \( Y_{it} = \beta_0 + \beta_1 X_{1,it} + \beta_2 X_{2,it} + \beta_3 X_{3,it} + \varepsilon_{it} \)
  • \( Y_{ro} = \beta_0 + \beta_1 X_{1,ro} + \beta_2 X_{2,ro} + \beta_3 X_{3,ro} + \varepsilon_{ro} \)
  • \( Y_{es} = \beta_0 + \beta_1 X_{1,es} + \beta_2 X_{2,es} + \beta_3 X_{3,es} + \varepsilon_{es} \)
  • etc

La modul și mai general, sunt folosiți doi indici, primul indexând unitatea statistică și al doilea predictorul:

  • \( Y_{1} = \beta_0 + \beta_1 X_{1,1} + \beta_2 X_{2,1} + \beta_3 X_{3,1} + \varepsilon_{1} \)
  • \( Y_{2} = \beta_0 + \beta_1 X_{1,2} + \beta_2 X_{2,2} + \beta_3 X_{3,2} + \varepsilon_{2} \)
  • \( Y_{3} = \beta_0 + \beta_1 X_{1,3} + \beta_2 X_{2,3} + \beta_3 X_{3,3} + \varepsilon_{3} \)
  • etc.

Dacă avem un termen liber, ca în acest caz, se poate standardiza scrierea lui β₀ pentru a fi un produs, ca și ceilalți termeni:

  • \( Y_{1} = \beta_0 \times 1 + \beta_1 \times X_{1,1} + \beta_2 \times X_{2,1} + \beta_3 \times X_{3,1} + \varepsilon_{1} \)
  • \( Y_{2} = \beta_0 \times 1 + \beta_1 \times X_{1,2} + \beta_2 \times X_{2,2} + \beta_3 \times X_{3,2} + \varepsilon_{2} \)
  • \( Y_{3} = \beta_0 \times 1 + \beta_1 \times X_{1,3} + \beta_2 \times X_{2,3} + \beta_3 \times X_{3,3} + \varepsilon_{3} \)
  • etc.

Se poate rescrie sistemul de ecuații de mai sus într-o singură formulă, numită ecuația matricială a modelului, \(\mathbf{Y} = \mathbf{X}\mathbf{\beta} + \mathbf{\epsilon}\).

În acea ecuație, \(\mathbf{X}\), numită matricea modelului (eng. design matrix), cuprinde toate valorile cunoscute ale predictorilor. Dacă modelul include un termen liber, i se va aloca și acestuia o primă coloană, cuprinzând numai valori 1:

Matricea \(\mathbf{Y}\) cuprinde toate valorile cunoscute ale variabilei previzionate, aici ale PIBului:

Matricile \(\mathbf{\beta}\) și \(\mathbf{\epsilon}\) vor cuprinde valorile parametrilor. Dacă dispunem doar de eșantioane, ecuația modelului devine \(\mathbf{Y} = \mathbf{X}\mathbf{\hat{\beta}} + \mathbf{\hat{\epsilon}}\), unde parametrii sunt înlocuiți de estimațiile lor:

  • \( \mathbf{\hat {\epsilon}} = \begin{pmatrix} \hat {\epsilon_1} \\ \hat {\epsilon_2} \\ \hat {\epsilon_3} \\ \hat {\epsilon_4} \\ \hat {\epsilon_5} \\ \hat {\epsilon_6} \end{pmatrix} \)
  • \( \mathbf{\hat {\beta}} = \begin{pmatrix} \hat {\beta_0} \\ \hat {\beta_1} \\ \hat {\beta_2} \\ \hat {\beta_3} \end{pmatrix} \)

Pentru ca suma reziduurilor să fie nulă, iar suma pătratelor reziduurilor să fie minimă, este necesar ca \( \hat{\mathbf{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y} \). Este posibil ca, la examene, să primiți direct valorile \( (\mathbf{X}^\top \mathbf{X})^{-1} \) și \( \mathbf{X}^\top \mathbf{y} \), situație în care va trebui să le înmulțiți, pentru a obține o matrice coloană cu β̂₀, β̂₁ etc.

Aici, \(\mathbf{X}^\top\) este transpusa lui X, iar \(\mathbf{X}^\top \mathbf{X}\), numită și matrice Gram a lui X, este o matrice pătrată, de dimensiuni k×k, unde k este numărul de predictori, incluzând termenul liber, dacă acesta este parte din model.

Inversul matricii Gram, \( (\mathbf{X}^\top \mathbf{X})^{-1}\), de aceleași dimensiuni ca matricea precedentă, este util în mai multe calcule:

În cursul unui examen scris, această matrice \( (\mathbf{X}^\top \mathbf{X})^{-1}\) este cel mai probabil să vă fie oferită în formă deja calculată.

Cealaltă matrice necesară estimării parametrilor este \( \mathbf{X}^\top \mathbf{Y}\), o matrice coloană cu k elemente, unde k este numărul de predictori, incluzând termenul liber, dacă acesta este parte din model:

Dacă aveți matricile \( (\mathbf{X}^\top \mathbf{X})^{-1}\) și \( \mathbf{X}^\top \mathbf{Y}\), veți putea calcula produsul, matricea-coloană cu toate cele k estimații ale parametrilor modelului, \(\mathbf{\hat{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{Y} \), unde k este numărul de predictori, incluzând termenul liber, dacă acesta este parte din model:

Se notează cu \(\mathbf{L}\) produsul (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \), pentru a putea scrie \(\mathbf{\hat{\beta}} = \mathbf{L} \mathbf{Y} \). Este mai puțin probabil ca la examen să trebuiască să calculați β̂ din L și Y, pentru că au dimensiuni mai mari (n coloane la prima și n linii la a doua, cu n – efectivul eșantionului).

Estimarea erorilor

Erorile sunt estimate de reziduuri. Din ecuația matricială a modelului reiese că reziduurile sunt \(\mathbf{\hat{\epsilon}} = \mathbf{Y} – \mathbf{X}\mathbf{\hat{\beta}}\).

Matricea \(\mathbf{X}\mathbf{\hat{\beta}}\) cuprinde, desigur, valorile estimate Ŷ, o matrice-coloană cu n elemente, unde n este efectivul eșantionului:

Deoarece \(\mathbf{\hat{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{Y} \), matricea valorilor estimate este \(\mathbf{\hat{Y}} = \mathbf{X} (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y} \).

Se notează \(\mathbf{H} = \mathbf{X} (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \), pentru ca să putem scrie \(\mathbf{\hat{Y}} = \mathbf{H} \mathbf{Y}\). Deși aceste trei matrici, Y, H și Ŷ, au dimensiuni mici, este puțin probabil ca testarea lor să fie subiect de examen. Totuși, relația Ŷ=HY este ușor de memorat.

Având matricile Ŷ și Y, putem calcula matricea reziduurilor, \(\mathbf{\hat{\epsilon}} = \mathbf{Y} – \mathbf{\hat{Y}}\), tot o matrice-coloană cu n elemente. Aici:

Estimarea SE pentru parametri

Abaterea medie standard (Standard Error, SE) este o măsură a incertitudinii pentru fiecare din estimații. În regresia liniară, estimațiile β̂₀, β̂₁ etc. vor avea, fiecare o abatere medie standard, cu care vom calcula intervale de încredere și vom evalua ipoteze statistice. Cuvântul “eroare” din sintagma “eroare standard” nu are legătură cu eroarea sau reziduul ce caracterizează fiecare observație (Y, X₁, X₂, X₃..), nefiind o măsură a imperfecțiunii modelului. În contextul matricilor:

  1. Suma pătratelor reziduurilor (SPR în română; SSresiduals în Excel) este o primă măsură a imperfecțiunilor modelului.
  2. Împărțirea lui SSresiduals la numărul său de grade de libertate (vezi capitolul despre intervale de încredere în regresia multiplă) estimează varianța (dispersia) erorilor, numită Mean Square Error (MSE) în Excel.
  3. Rădăcina pătrată a MSE este abaterea standard a regresiei (standard error of the estimate), o măsură versatilă a incertitudinii modelului, de asemenea detaliată in capitolul despre intervale de încredere în regresia multiplă.

În acest exemplu, tabelul ANOVA și abaterea standard a regresiei sunt:

La modul general, abaterile pătratice standard necesare estimării parametrilor ca intervale de încredere și testării de ipoteze se obțin din matricea \( \hat \sigma^2 (\mathbf{X}^\top\mathbf{X})^{-1}\), unde X este matricea modelului. Această matrice, cu (p+1) x (p+1) elemente, unde p este numărul de predictor excluzând termenul liber, cuprinde pe diagonală varianțele fiecărei estimări:

  • al doilea element al diagonalei, notat c₁₁, este Var[β̂₁]
  • al treilea element al diagonalei, notat c₂₂, este Var[β̂₂]
  • etc.

Celelalte valori din matricea \( \hat \sigma^2 (\mathbf{X}^\top\mathbf{X})^{-1}\) sunt covarianțe ale estimărilor, ca de exemplu cov(β̂₁, β̂₂) etc. Din acest motiv, matricea \( \hat \sigma^2 (\mathbf{X}^\top\mathbf{X})^{-1}\) se mai numește matricea de varianță-covarianță, si se notează \(Var[ \mathbf{\hat \beta}]\). Aici, folosind \((\mathbf{X}^\top\mathbf{X})^{-1}\) de mai sus, matricea de varianță-covarianță este:

Rădăcina pătrată a numerelor de pe diagonala lui \(Var[ \mathbf{\hat \beta}]\) ne va reda abaterea standard (SE) pentru fiecare din estimatorii parametrilor modelului, SE[β̂₁], SE[β̂₂] etc., Atenție: la unele teste, vi se va da \( \hat \sigma^2 (\mathbf{X}^\top\mathbf{X})^{-1}\) și la altele \( (\mathbf{X}^\top\mathbf{X}) ^ {-1}\), în al doilea caz fiind necesar să înmulțiți elementele de pe diagonală cu MSE=σ̂² pentru a obține erorile standard. Similar, este posibil să primiți MSE sau σ̂.

Aici, șirul SE este:

Un caz mai complicat este cel al termenului liber, pentru care în loc de primul element de pe diagonală, se folosește

\(c_{00} \;=\; \frac{1}{n} \;+\; \bar{\mathbf{x}}^{\top}\,\bigl(X_X^{\top}X_X\bigr)^{-1}\,\bar{\mathbf{x}} \)

unde \( \bar{\mathbf{x}} \) este matricea cu mediile fiecărei coloane din \( \mathbf{X} \), iar Xₓ este matricea modelului fără coloana cu valori 1 corespunzătoare termenului liber.

Previzionare din regresia liniară cu calcul matricial

Dat fiind setul de valori X₀₁, X₀₂, X₀₃ etc pentru o unitate statistică nouă, se poate construi o matrice \( \mathbf{X_0} \), cu o singură linie din acestea, la care trebuie atașată 1 pe coloana inițială, dacă modelul include termen liber. În acest caz

  • estimarea punctuală este \(\hat{Y}_0 = \hat{\beta} X_0\)
  • pentru medii, varianța estimației este \(\mathrm{Var}(\widehat{\overline{Y}}_0) = \hat{\sigma}^2\, X_0^{\top} (X^{\top} X)^{-1} X_0\)
  • pentru valori individuale, varianța estimației este \(\mathrm{Var}(\hat{Y}_0) = \hat{\sigma}^2\, (1 + X_0^{\top} (X^{\top} X)^{-1} X_0)\).

Recapitulare

  • În regresiei liniare, ecuațiile din modelul populațional se pot scrie sub formă de matrici, sub forma \(\mathbf{Y} = \mathbf{X}\mathbf{\beta} + \mathbf{\epsilon}\), unde \(\mathbf{Y}\), \(\mathbf{\beta}\) și \(\mathbf{\epsilon}\) sunt matrici coloană cu toate valorile variabileo previzionate, ale parametrilor și respectiv ale erorilor. În cazul eșantioanelor, matricile de parametri sunt înlocuite de matrici cu estimații, ecuația modelului devenind \(\mathbf{Y} = \mathbf{X}\mathbf{\hat{\beta}} + \mathbf{\hat{\epsilon}}).
  • Matricea \(\mathbf{X}\), numită matricea modelului, cuprinde o linie pentru fiecare unitate statistică (observație) și o coloană pentru fiecare predictor. Dacă modelul include un termen liber, coloana predictorului este prima și cuprinde numai valori 1.
  • Metoda celor mai mici pătrate necesită ca suma pătratelor erorilor să fie minimizată, ceea ce se obține când \( \hat{\mathbf{\beta}}\) este estimat cu \( (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y} \).
  • Matricea valorilor ajustate \( \hat{\mathbf{Y}}\) va fi \(\mathbf{X}\mathbf{\hat{\beta}}\), iar matricea reziduurilor \(\mathbf{\epsilon}\) va fi \( \mathbf{Y} – \hat{\mathbf{Y}}\).
  • Suma pătratelor reziduurilor, MSE (raportul SSresiduals / dfresiduals) și abaterea standard a regresiei (σ̂, adică √MSE) sunt măsuri ale imperfecțiunii modelului, ultimele două sintetizând dimensiunea erorilor (nu a reziduurilor!).
  • Matricea \( \hat \sigma^2 (\mathbf{X}^\top\mathbf{X})^{-1}\), numită matrice varianță-covarianță, cuprinde pe diagonala sa valori necesare calculării abaterii standard a fiecărei estimații de parametru. Dacă modelul include un termen liber, primul număr de pe diagonală necesită calcule suplimentare relativ prelungite, însă, în restul cazurilor, numerele de pe diagonală sunt varianțele estimațiilor. În cazul includerii termenului liber în model, al doilea număr de pe diagonală este Var[β₁], adică pătratul lui SE[β₁]; al treilea număr de pe diagonală este Var[β₂] = (SE[β₂])² etc.

Problemă propusă

Cunoaștem următoarele informații despre un eșantion de state:

Să se estimeze parametrii unui model de regresie liniară multiplă PIB = β₀ + β₁×Pop + β₂ × RPFM, precum și abaterile standard ale acestora.

Răspuns