La linea de mínimos cuadrados

El método más usado para estimar el modelo lineal (la línea que mejor describe los puntos ajustados) es llamado mínimos cuadrados, para este punto recordemos que la ecuación de la línea recta es: \(y = m x + b_0\).

Expresiones para estimar la línea recta de mínimos cuadrados

Las ecuaciones para estimar la línea recta por mínimos cuadrados ara un conunto de datos \(Y = y_1, y_2, \dots, y_n\) y \(X = x_1, x_2, \dots, x_n\), con promedios \(\bar{y}\) y \(\bar{x}\), respectivamente y \(\hat{\beta}_0\) y \(\hat{\beta}_1\) son los valores estimados para la intersección con el eje \(y\), \(\beta_0\) y la pendiente de la recta \(\beta_1\).

\[ \hat{\beta}_1 = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n (x_i - \bar{x})^2} \]

\[ \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x} \]

Ejemplo Deprès et al. establecieron que la topografía del tejido adiposo (AT) está asociada con complicaciones metabólicas consideradas como un riesgo par la enfermedad cardiovascular. Es importante, establecieron en dicho trabajo, medir la cantidad de tejido adiposo itraabdominal como parte de la evaluación de riesgo individual de enfermedad cardiovascular. Tomografía computarizada (CT), la única técnica que puede medir de forma precisa y confiable la cantidad de AT intraabdominal profunda, sin embargo, es costosa y expone al sujeto a radiación ionizante (rayos X), además es una técnica que no es accesible a todos los médicos. Deprès y sus colegas condujeron su estudio para desarrollar ecuaciones que predigan la cantidad de AT intraabdominal profunda a partir de medidas antropométricas simples. Sus sujetos de estudio fueron hombres adultos de entre 18 y 42 años de edad, sin enfermedad de desorden metabólico que requiriera medicación. Entre las medidas tomadas a cada sujeto estuvieron la de TA intraabdominal profunda por CT y la circunferencia de la cintura, listadas en el archivo EXA_C09_S03_01.csv. La pregunta es; que tan bien se puede estimar la cantidad de TA intraabdominal profunda partir se la medida de la circunferencia de la cintura. Como la medida de AT intraabdominal profunda es la variable a evaluar y hacer predicciones, se considera la variable dependiente. La variable de la circunferencia de la cintura se tomará como la medida independiente.

Waist <- read_csv("~/Dropbox/GitHub/ProbEstad/DataSets/ch09_all/EXA_C09_S03_01.csv", show_col_types = FALSE)
Waist

# A tibble: 109 × 3
    SUBJ     X     Y
   <dbl> <dbl> <dbl>
 1     1  74.8  25.7
 2     2  72.6  25.9
 3     3  81.8  42.6
 4     4  84.0  42.8
 5     5  74.6  29.8
 6     6  71.8  21.7
 7     7  80.9  29.1
 8     8  83.4  33.0
 9     9  63.5  11.4
10    10  73.2  32.2
# … with 99 more rows

plot(Y ~ X, data = Waist, pch = 20, xlab = "Circunferencia de cintura (cm)", ylab = "área de TA intraabdominal profunda")

# mean lines for y and x.
abline(v = mean(Waist$X), col = 2, lty = 2)
abline(h = mean(Waist$Y), col = 2, lty = 2)
# Estimation of the linear model
Waist_lm <- lm(Y ~ X, data = Waist)
summary(Waist_lm)

Call:
lm(formula = Y ~ X, data = Waist)

Residuals:
     Min       1Q   Median       3Q      Max 
-107.288  -19.143   -2.939   16.376   90.342 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -215.9815    21.7963  -9.909   <2e-16 ***
X              3.4589     0.2347  14.740   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 33.06 on 107 degrees of freedom
Multiple R-squared:   0.67, Adjusted R-squared:  0.667 
F-statistic: 217.3 on 1 and 107 DF,  p-value: < 2.2e-16

# plot of the estimated linear model
abline(a = -215.98, b = 3.4589, col = 6)

Estimación del Coeficiente de Determinación de la Población

Siempre es importante preguntarnos si la recta ajustada es una buena representación de los datos, es decir que tan buen trabajo hizo el ajuste al modelo lineal. Para es to evaluamos la ecuación de la regresión basada en los datos de la muestra (\(\tilde{y}_i\)), el resultado es llamado el coeficiente de determinación de la muestra, \(r^2\).

\[ r^2 = \frac{\sum(\hat{y}_i - \bar{y})^2}{\sum(y_i - \bar{y})^2} \]

Esta R-cuadrada es calculada directamente de los cuadrados de la diferencias de la \(y_i\) estimadas por el modelo en cada punto \(x_i\) comparados con los cuadrados de los datos originales \(y_i\) (en R es llamada Miltiple R-squared). El coeficiente de determinación de la muestra mide la “calidad” (que tan cerca) está la ecuación de regresión a los valores observados originales.

El coeficiente de determinación de la muestra es una estimado puntual de \(\rho^2\) el Coeficiente de Determinación de la Población. El Coeficiente de determinación de la población \(\rho^2\) tiene la misma función relativa a la población que \(r^2\) tiene a la muestra. Este coeficiente, muestra cual proporción de la variación de la población en \(Y\) es explicada por la regresión de \(Y\) en \(X\). Un estimador sin sesgo de \(\rho^2\) está dado por

\[ \tilde{r}^2 = 1 - \frac{\sum(y_i - \hat{y}_i)^2/(n-2)}{\sum(y_i - \bar{y})^2/(n-1)} \] Esta es la cantidad llamada R-sq (en R es la llamada Adjusted R-squared).

EjemploLos escores listados en EXR_C09_S03_02.csv, muestran las evaluaciones de la enfermera (X) y del médico (Y) para la condición de 10 pacientes al momento de admisión de un centro de trauma. Muestra una gráfica de los puntos y estime la regresión lineal.

Solución primero resolveremos el modelo en forma general estimando los dos parámetros de ajuste de mínimos cuadrados, el término independiente (intersección con el eje \(y\)) \(\beta_0\) y el segundo término o pendiente \(\beta_1\) con el modelo lineal Y ~ X. Posteriormente estimaremos el modelo “sin el término independiente”, forzaremos que el modelo en R estime sin el coeficiente \(\beta_0\), Y ~ X - 1. En la gráfica el modelo lineal completo es la línea roja y el modelo que pasa por el origen es la línea azul.

Scores <- read_csv("~/Dropbox/GitHub/ProbEstad/DataSets/ch09_all/EXR_C09_S03_02.csv", show_col_types = FALSE)
Scores

# A tibble: 10 × 2
       X     Y
   <dbl> <dbl>
 1    18    23
 2    13    20
 3    18    18
 4    15    16
 5    10    14
 6    12    11
 7     8    10
 8     4     7
 9     7     6
10     3     4

# X enfermera, Y médico
plot(Y ~ X, data = Scores, pch = 20, xlab = "Eval. enfermera", ylab = "Eval. médico", xlim = c(0, 20), ylim = c(0, 25))
abline(v = mean(Scores$X), col = 2, lty = 2)
abline(h = mean(Scores$Y), col = 2, lty = 2)

# MinSqr 
Scores_lm <- lm(Y ~ X, data = Scores)
summary(Scores_lm)

Call:
lm(formula = Y ~ X, data = Scores)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1988 -2.3808 -0.1638  1.8393  4.7189 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.2112     2.0557   0.589 0.571993    
X             1.0823     0.1723   6.283 0.000237 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.765 on 8 degrees of freedom
Multiple R-squared:  0.8315,    Adjusted R-squared:  0.8104 
F-statistic: 39.47 on 1 and 8 DF,  p-value: 0.0002372

# Good example for forcing the model to set beta_0 =0
Scores_lm_0 <- lm(Y ~ X - 1, data = Scores)
summary(Scores_lm_0)

Call:
lm(formula = Y ~ X - 1, data = Scores)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1348 -2.0674  0.5421  2.1601  4.7360 

Coefficients:
  Estimate Std. Error t value Pr(>|t|)    
X  1.17416    0.07056   16.64 4.57e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.663 on 9 degrees of freedom
Multiple R-squared:  0.9685,    Adjusted R-squared:  0.965 
F-statistic: 276.9 on 1 and 9 DF,  p-value: 4.566e-08

# plot of the linear model for Scores
abline(a = Scores_lm$coefficients[1], b = Scores_lm$coefficients[2], col = 2)
abline(a = 0, b = Scores_lm_0$coefficients, col = 4)

# Error plots for lm in R
plot(Scores_lm)

# Error plots for the second model
plot(Scores_lm_0)

Suposiciónes del modelo de correlación

A continuación listamos las suposiciones necesarias para hacer inferencias de una población a partir de muestreo de una distribución bivariada:
1. Para cada valor de \(X\) hay una subpoblación normalmente distribuida de valores de \(Y\).
2. Para cada valor de \(Y\) hay una subpoblación normalmente distribuida de valores de \(X\).
3. La distribución conjunta de \(X\) y \(Y\) es una distribución normal bivariada.
4. Las subpoblaciones de valores de \(Y\) tienen la misma varianza.
5. Las subpoblaciones de valores de \(X\) tienen la misma varianza.

El coeficiente de correlación

La distribución bivariada normal que describe estos datos tiene cinco parámetros: \(\sigma_x, \sigma_y, \mu_x, \mu_y,\) y \(\rho\). Los cuatro primeros son las varianza de \(X\) y \(Y\) y sus medias respectivas y \(\rho\) es el coeficiente de correlación de la población, que mide la fuerza de la relación lineal entre las variables \(Y\) y \(X\). Por lo general estamos interesados en determinar si \(\rho \neq 0\), esto es. si \(Y\) y \(X\) están linealmente correlacionadas. Como \(\rho\) es por lo general desconocido, es necesario extraer muestras de las poblaciones para estimas el valor de \(r\) y determinar el valor de \(\rho\). Una hipótesis nula más utilizada es \(H_0: \rho = 0 \space | H_A: \rho \neq 0\).
El valor de r se puede estimar por:

\[ r = \frac{n \sum x_i y_i - (\sum x_i)(\sum y_i)}{\sqrt{n \sum x_i^2 - (\sum x_i)^2} \sqrt{n \sum y_i^2 - ( \sum y_i)^2}} \]

Ejemplo El propósito del estudio de Kwast-Rabben y colaboradores, fue analizar los potenciales somatosensoriales evocados (SEPs) y su interrelación seguida de la estimulación de los dígitos I, III y V de la mano. Los investigadores querían establecer un antecedente para criterio de referencia en una población control. Por lo tanto se reclutaron sujetos control sanos para el estudio. En un futuro este estudio podría ser importante ya que los SEPs control se podrían usar para establecer disturbios funcionales en pacientes con sospecha de lesión en la raíz cervical que presentan dolor, y síntomas sensoriales. En el estudio se estimularon los dedos con una intensidad por debajo del punto de dolor. Las respuestas medulares se registraron con electrodos fijos adheridos con crema adherente a la piel del sujeto. Una de las relaciones de interés fue la correlación entre la estatura del sujeto (cm) y la latencia medular pico (Cv) del SEP. Los datos obtenidos se muestran en EXA_C09_S07_01.csv.

SEP <- read_csv("~/Dropbox/GitHub/ProbEstad/DataSets/ch09_all/EXA_C09_S07_01.csv", show_col_types = FALSE)
plot(CV ~ HEIGHT, data = SEP, pch = 20, xlab = "Estatura (cm)", ylab = "Latencia pico (Cv)")

# cor.test(HEIGHT, CV, alternative = c("two.sided"), method = c("pearson"), conf.level = 0.95)
cor.test(~ HEIGHT + CV, data = SEP, alternative = c("two.sided"), method = c("pearson"), conf.level = 0.95)

    Pearson's product-moment correlation

data:  HEIGHT and CV
t = 19.781, df = 153, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7967316 0.8869721
sample estimates:
      cor 
0.8478829 

El propósito del estudio de Brown y Presley fue caracterizar la hepatitis A aguda, en pacientes mayores de 40 años. Revisaron en forma retrospectiva los registros de 20 pacientes que fueron diagnosticados con hepatitis A aguda, que no fueron hospitalizados. De interés fue el hecho que se usó la edad (en años) para predecir los niveles de bilirrubina (mg/dl). Los datos colectados están en EXR_C09_S07_01.csv.

Hepa <- read_csv("~/Dropbox/GitHub/ProbEstad/DataSets/ch09_all/EXR_C09_S07_01.csv", show_col_types = FALSE)
plot(BILI ~ AGE, data = Hepa, pch = 20, xlab = "Edad (años)", ylab = "Nivel de Bilirrubina (mg/dl)")
cor.test(~ AGE + BILI, data = Hepa, alternative = c("two.sided"), method = c("pearson"), conf.level = 0.95)

    Pearson's product-moment correlation

data:  AGE and BILI
t = 2.2331, df = 18, p-value = 0.03848
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.02927797 0.75306959
sample estimates:
      cor 
0.4657643 

# Linear model
Hepa_lm <- lm(BILI ~ AGE, data = Hepa)
abline(a = Hepa_lm$coefficients[1], b = Hepa_lm$coefficients[2], col = 2)

summary(Hepa_lm)

Call:
lm(formula = BILI ~ AGE, data = Hepa)

Residuals:
   Min     1Q Median     3Q    Max 
-7.203 -2.927 -1.529  2.812  9.263 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -2.9842     4.5962  -0.649   0.5244  
AGE           0.1793     0.0803   2.233   0.0385 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.414 on 18 degrees of freedom
Multiple R-squared:  0.2169,    Adjusted R-squared:  0.1734 
F-statistic: 4.987 on 1 and 18 DF,  p-value: 0.03848