MedCalcEjemplo de script: Regresión múltiple mediante aritmética matricial
Este ejemplo muestra la formulación matricial de regresión múltiple e ilustra cómo calcularla utilizando el script MedCalc.
En general, la formulación matricial de la regresión múltiple es $ Y = Xb + e $ donde Y es un vector que contiene los datos dependientes, X es una matriz que contiene los datos x independientes y e es un vector que contiene los residuos.
En el ejemplo, utilizaremos los siguientes datos de muestra:
Script de MedCalc:
<data> varx1 varx2 vary 2.3 2.0 11.2 5.1 2.1 18.8 6.9 4.3 31.6 9.0 5.8 4.2 12.0 3.2 5.1 18.3 6.2 1.1 21.6 8.0 1.6 </data>
Vector de variable dependiente Y
Y es una matriz (n,1) que contiene los valores y de la variable dependiente:
$$ Y=\begin{bmatrix}
y_1 \\
y_2 \\
\vdots \\
y_{n} \end{bmatrix}$$
Script de MedCalc:
y={vary}; // Creates a matrix y from the data variable vary
Matriz aumentada X
X es una matriz (n,k+1) que contiene una columna de 1 y columnas de las distintas variables x:
$$ X=\begin{bmatrix}1 & x_{11}&x_{12}&\vdots&x_{1k}\\
1 & x_{21}&x_{22}&\vdots&x_{2k}\\
1 & \vdots&\vdots&\vdots&\vdots\\
1 & x_{n1}&x_{n2}&\vdots&x_{nk}\\
\end{bmatrix} $$
Script de MedCalc:
x={varx1,varx2}; // Creates a matrix x from the data variables varx1 and varx2
// counts
n = nrow(x);
k = ncol(x);
// add a column with 1s for the constant
// gives the augmented matrix xa
c=ones(n,1);
xa=cbind(c,x);
Vector de coeficientes de regresión b
b es una matriz (k+1,1) que contiene los coeficientes de regresión
$$ b = \begin{bmatrix}b_0 \\
b_1 \\
\vdots \\
b_{k} \end{bmatrix} $$
Vector de residuos e
e es una matriz (n,1) que contiene los residuos
$$ e=\begin{bmatrix}e_1 \\
e_2 \\
\vdots \\
e_{n} \end{bmatrix}$$
Computación b
$$ b = (X^{'}X)^{-1}X^{'}Y $$Script de MedCalc:
xt=transpose(xa); xtx=xt*xa; xtxI=inv(xtx); b=xtxI*xt*y;
Cálculo de valores predichos
La matriz (n,1) $ \hat{Y} $ contiene los valores predichos:
$$ \hat{Y} = Xb $$Script de MedCalc:
ypred=xa*b;
Cálculo de los residuos e
$$ e = Y - \hat{Y} $$Script de MedCalc:
e=y-ypred;
Cálculo del coeficiente de determinación R 2
$$ R^2 = y^{'}Xb(y^{'}y)^{-1} $$usando
$$ y = Y - \bar{Y} $$Script de MedCalc:
ya=y-avg(y); yat=transpose(ya); r=yat*xa*b*inv(yat*y);
Suma total de cuadrados
$$ Total = y^{'}y $$Script de MedCalc:
sstotal=yat*y;
Sumas residuales de cuadrados
$$ SSresid = e^{'}e $$Script de MedCalc:
et=transpose(e); ssresid=et*e;
Regresión suma de cuadrados
$$ SSregres = SStotal - SSresid $$Script de MedCalc:
ssregres=sstotal-ssresid;
Covarianza y errores estándar de los coeficientes de regresión
$$ C = \frac{SSresid}{nk-1} (X^{'}X)^{-1} $$Las raíces cuadradas de las diagonales de C son los errores estándar de los coeficientes de regresión.
Script de MedCalc:
C=ssresid/(n-k-1); C=C*inv(transpose(xa)*xa); C=sqrt(C); se=diag(C);
Calcular valores T y valores F
Script de MedCalc:
T=b&/se; // elementwise division pT=TDIST(T,n-k-1); F=(ssregres/k)/(ssresid/(n-k-1)); pF=FDIST(f,k,n-k-1);
Residuos de prueba para distribución normal
Script de MedCalc:
pR=testnorm(shapwilk,e);
Crear informe
Script de MedCalc:
print 'Coefficient of determination R²: ',r:4; println ; println ; // print 'Coefficients, Standard errors, T, P ' ; println ; print 'constant: ',b[1]:4:8,se[1]:4:8,T[1]:4:8,pT[1]:4:8; println ; print 'b1: ',b[2]:4:8,se[2]:4:8,T[2]:4:8,pT[2]:4:8; println ; print 'b2: ',b[3]:4:8,se[3]:4:8,T[3]:4:8,pT[3]:4:8; println ; println ; // print 'Analysis of variance' ; println ; print 'SStotal: ',sstotal:4:8; println ; print 'SSregres: ',ssregres:4:8; println ; print 'SSresid: ',ssresid:4:8; println ; println ; // print 'F-ratio: ',F:4:8; println ; print 'P : ',pF:4:8; println ; println ; // print 'Test residuals for normal distribution P = ',pR:4; println ;
La salida del script
Este script crea la siguiente salida:
Coefficient of determination R²: 0.4069262 Coefficients, Standard errors, T, P constant: 19.6395 9.6496 2.0353 0.1115 b1: -1.2889 1.2139 -1.0618 0.3482 b2: 1.0459 3.7869 0.2762 0.7961 Analysis of variance SStotal: 751.0086 SSregres: 305.6051 SSresid: 445.4035 F-ratio: 1.3723 P : 0.3517 Test residuals for normal distribution P = 0.4909
El script completo
// regional settings
$decimalsymbol= '.' ;
$listseparator= ',' ;
//
cls ;
clear data;
clear mem;
<data>
varx1 varx2 vary
2.3 2.0 11.2
5.1 2.1 18.8
6.9 4.3 31.6
9.0 5.8 4.2
12.0 3.2 5.1
18.3 6.2 1.1
21.6 8.0 1.6
</data>
// create matrices from data
y={vary};
x={varx1,varx2};
// counts
n = nrow(x);
k = ncol(x);
// add a column with 1s for the constant
// gives the augmented matrix xa
c=ones(n,1);
xa=cbind(c,x);
// b = (X’X)^-1 X’ Y
xt=transpose(xa);
xtx=xt*xa;
xtxI=inv(xtx);
// calculate coefficients b
b=xtxI*xt*y;
// calculate predicted values
ypred=xa*b;
// calculate residuals
e=y-ypred;
//Coefficient of determination R²
ya=y-avg(y);
yat=transpose(ya);
r=yat*xa*b*inv(yat*y);
// Total sum of squares
sstotal=yat*y;
// Residual sum of squares;
et=transpose(e);
ssresid=et*e;
// Regression sum of squares
ssregres=sstotal-ssresid;
// Covariance and
// Standard Errors of regression coefficients
C=ssresid/(n-k-1);
C=C*inv(transpose(xa)*xa);
C=sqrt(C);
se=diag(C);
// T-values and F-value;
T=b&/se; // elementwise division
pT=TDIST(T,n-k-1);
F=(ssregres/k)/(ssresid/(n-k-1));
pF=FDIST(f,k,n-k-1);
// Test residuals for normal distribution
pR=testnorm(shapwilk,e);
// report
print 'Coefficient of determination R²: ',r:4; println ;
println ;
//
print 'Coefficients, Standard errors, T, P ' ; println ;
print 'constant: ',b[1]:4:8,se[1]:4:8,T[1]:4:8,pT[1]:4:8; println ;
print 'b1: ',b[2]:4:8,se[2]:4:8,T[2]:4:8,pT[2]:4:8; println ;
print 'b2: ',b[3]:4:8,se[3]:4:8,T[3]:4:8,pT[3]:4:8; println ;
println ;
//
print 'Analysis of variance' ; println ;
print 'SStotal : ',sstotal:4:8; println ;
print 'SSregres: ',ssregres:4:8; println ;
print 'SSresid : ',ssresid:4:8; println ;
println ;
//
print 'F-ratio: ',F:4:8; println ;
print 'P : ',pF:4:8; println ;
println ;
//
print 'Test residuals for normal distribution P = ',pR:4; println ;
Digimizer