1 Introduccion

El presente documento desarrolla un analisis estadistico completo sobre el conjunto de datos Salaries, disponible en el paquete car de R. Este conjunto de datos contiene informacion sobre los salarios de profesores universitarios en Estados Unidos durante el ano academico 2008-2009, y es ampliamente utilizado en la ensenanza de modelos de regresion lineal multiple.

El analisis esta estructurado de la siguiente manera:

  1. Descripcion y exploracion del conjunto de datos
  2. Analisis exploratorio de datos (EDA) con visualizaciones
  3. Evaluacion de supuestos: normalidad y correlacion
  4. Analisis de multicolinealidad previo al modelo
  5. Seleccion de variables con el metodo Forward (funcion step)
  6. Ajuste e interpretacion del modelo de regresion lineal multiple
  7. Analisis completo de residuos
  8. Deteccion de valores atipicos e influyentes
  9. Analisis de multicolinealidad del modelo final

2 Carga de paquetes y datos

# Instalacion de paquetes si no estan disponibles
paquetes <- c("car", "ggplot2", "dplyr", "nortest", "lmtest",
              "corrplot", "gridExtra", "knitr", "kableExtra")

for (p in paquetes) {
  if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
}

library(car)
library(ggplot2)
library(dplyr)
library(nortest)
library(lmtest)
library(corrplot)
library(gridExtra)
library(knitr)
library(kableExtra)

# Carga del conjunto de datos
data("Salaries", package = "car")

# Vista rapida
head(Salaries, 10) |>
  kable(caption = "Primeras 10 observaciones del conjunto de datos Salaries") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Primeras 10 observaciones del conjunto de datos Salaries
rank discipline yrs.since.phd yrs.service sex salary
Prof B 19 18 Male 139750
Prof B 20 16 Male 173200
AsstProf B 4 3 Male 79750
Prof B 45 39 Male 115000
Prof B 40 41 Male 141500
AssocProf B 6 6 Male 97000
Prof B 30 23 Male 175000
Prof B 45 45 Male 147765
Prof B 21 20 Male 119250
Prof B 18 18 Female 129000

3 Descripcion del Conjunto de Datos

3.1 Estructura general

str(Salaries)
## 'data.frame':    397 obs. of  6 variables:
##  $ rank         : Factor w/ 3 levels "AsstProf","AssocProf",..: 3 3 1 3 3 2 3 3 3 3 ...
##  $ discipline   : Factor w/ 2 levels "A","B": 2 2 2 2 2 2 2 2 2 2 ...
##  $ yrs.since.phd: int  19 20 4 45 40 6 30 45 21 18 ...
##  $ yrs.service  : int  18 16 3 39 41 6 23 45 20 18 ...
##  $ sex          : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 1 ...
##  $ salary       : int  139750 173200 79750 115000 141500 97000 175000 147765 119250 129000 ...

3.2 Descripcion de las variables

El conjunto de datos contiene 397 observaciones y 6 variables. A continuacion se describe cada variable:

Variable Tipo Descripcion
rank Categorica Rango academico del profesor: AsstProf (Profesor Asistente), AssocProf (Profesor Asociado), Prof (Profesor Titular)
discipline Categorica Area disciplinar: A (Disciplinas teoricas, ej. humanidades) o B (Disciplinas aplicadas, ej. ciencias)
yrs.since.phd Numerica Anos transcurridos desde la obtencion del doctorado
yrs.service Numerica Anos de servicio en la institucion
sex Categorica Sexo del profesor: Female o Male
salary Numerica Salario anual en dolares (variable respuesta)

Variable respuesta: salary (salario anual en dolares americanos).

Variables explicativas: rank, discipline, yrs.since.phd, yrs.service, sex.

3.3 Rangos y categorias de las variables categoricas

Es fundamental entender los niveles de las variables categoricas para interpretar correctamente los resultados del modelo:

# rank
cat("Variable: rank\n")
## Variable: rank
table(Salaries$rank)
## 
##  AsstProf AssocProf      Prof 
##        67        64       266
cat("\nVariable: discipline\n")
## 
## Variable: discipline
table(Salaries$discipline)
## 
##   A   B 
## 181 216
cat("\nVariable: sex\n")
## 
## Variable: sex
table(Salaries$sex)
## 
## Female   Male 
##     39    358

Descripcion detallada de los niveles:

  • rank: Jerarquia academica ascendente
    • AsstProf: Profesor Asistente (nivel inicial de la carrera academica, n = 67)
    • AssocProf: Profesor Asociado (nivel intermedio, n = 64)
    • Prof: Profesor Titular (nivel maximo, n = 266)
  • discipline: Tipo de disciplina academica
    • A: Disciplinas teoricas o humanisticas (n = 181)
    • B: Disciplinas aplicadas o cientificas (n = 216)
  • sex: Sexo del docente
    • Female: Mujer (n = 39)
    • Male: Hombre (n = 358)

4 Analisis Exploratorio de Datos

4.1 Estadisticas descriptivas

summary(Salaries) |> print()
##         rank     discipline yrs.since.phd    yrs.service        sex     
##  AsstProf : 67   A:181      Min.   : 1.00   Min.   : 0.00   Female: 39  
##  AssocProf: 64   B:216      1st Qu.:12.00   1st Qu.: 7.00   Male  :358  
##  Prof     :266              Median :21.00   Median :16.00               
##                             Mean   :22.31   Mean   :17.61               
##                             3rd Qu.:32.00   3rd Qu.:27.00               
##                             Max.   :56.00   Max.   :60.00               
##      salary      
##  Min.   : 57800  
##  1st Qu.: 91000  
##  Median :107300  
##  Mean   :113706  
##  3rd Qu.:134185  
##  Max.   :231545

Las estadisticas descriptivas revelan lo siguiente:

  • El salario promedio es de 113,706 dolares anuales, con un minimo de 57,800 y un maximo de 231,545.
  • Los anos desde el doctorado van de 1 a 56 anos, con una media de 22.3.
  • Los anos de servicio van de 0 a 60 anos.
  • La muestra esta compuesta mayoritariamente por hombres (90.2%) y por Profesores Titulares (67%).

4.2 Distribucion de la variable respuesta: salary

par(mfrow = c(1, 2))

hist(Salaries$salary,
     main   = "Histograma del Salario",
     xlab   = "Salario (USD)",
     ylab   = "Frecuencia",
     col    = "steelblue",
     border = "white",
     breaks = 25)
abline(v = mean(Salaries$salary), col = "red", lwd = 2, lty = 2)
legend("topright", legend = "Media", col = "red", lty = 2, lwd = 2)

boxplot(Salaries$salary,
        main   = "Boxplot del Salario",
        ylab   = "Salario (USD)",
        col    = "steelblue",
        border = "gray30")

par(mfrow = c(1, 1))

El histograma muestra una distribucion con sesgo hacia la derecha (asimetria positiva), lo cual es comun en datos de salarios. Hay presencia de valores elevados que desplazan la media por encima de la mediana. El boxplot confirma la existencia de valores atipicos en la parte superior de la distribucion.

4.3 Distribucion de variables numericas continuas

par(mfrow = c(1, 2))

hist(Salaries$yrs.since.phd,
     main   = "Anos desde el Doctorado",
     xlab   = "Anos",
     ylab   = "Frecuencia",
     col    = "darkorange",
     border = "white",
     breaks = 20)

hist(Salaries$yrs.service,
     main   = "Anos de Servicio",
     xlab   = "Anos",
     ylab   = "Frecuencia",
     col    = "darkgreen",
     border = "white",
     breaks = 20)

par(mfrow = c(1, 1))

Ambas variables numericas presentan distribuciones asimetricas. yrs.since.phd muestra mayor dispersion, mientras que yrs.service tiene una concentracion de valores en los primeros anos.

4.4 Boxplots de Salario segun variables categoricas

Esta seccion es fundamental para el analisis exploratorio: permite visualizar como se distribuye el salario dentro de cada categoria de las variables explicativas, e identificar diferencias entre grupos antes de ajustar el modelo.

4.4.1 Salario segun Rango Academico (rank)

# Orden logico de menor a mayor jerarquia
Salaries$rank <- factor(Salaries$rank,
                         levels = c("AsstProf", "AssocProf", "Prof"))

ggplot(Salaries, aes(x = rank, y = salary, fill = rank)) +
  geom_boxplot(outlier.colour = "red", outlier.shape = 16, outlier.size = 2) +
  stat_summary(fun = mean, geom = "point", shape = 23,
               size = 3, fill = "white", color = "black") +
  scale_fill_manual(values = c("AsstProf"  = "#4C72B0",
                                "AssocProf" = "#DD8452",
                                "Prof"      = "#55A868")) +
  scale_x_discrete(labels = c("AsstProf"  = "Profesor\nAsistente",
                               "AssocProf" = "Profesor\nAsociado",
                               "Prof"      = "Profesor\nTitular")) +
  scale_y_continuous(labels = scales::comma) +
  labs(title  = "Distribucion del Salario segun Rango Academico",
       subtitle = "El rombo blanco indica la media de cada grupo",
       x = "Rango Academico",
       y = "Salario Anual (USD)") +
  theme_bw() +
  theme(legend.position = "none",
        plot.title    = element_text(face = "bold", size = 13),
        plot.subtitle = element_text(size = 10, color = "gray40"))

Interpretacion: Se observa una tendencia clara y esperada: el salario aumenta con el rango academico. Los Profesores Titulares (Prof) tienen la mayor mediana salarial y la mayor dispersion, lo que indica heterogeneidad dentro de ese grupo. Los Profesores Asistentes (AsstProf) presentan la menor dispersion y salarios mas homogeneos. Los puntos rojos corresponden a valores atipicos dentro de cada categoria.

4.4.2 Salario segun Disciplina (discipline)

ggplot(Salaries, aes(x = discipline, y = salary, fill = discipline)) +
  geom_boxplot(outlier.colour = "red", outlier.shape = 16, outlier.size = 2) +
  stat_summary(fun = mean, geom = "point", shape = 23,
               size = 3, fill = "white", color = "black") +
  scale_fill_manual(values = c("A" = "#C44E52", "B" = "#8172B2")) +
  scale_x_discrete(labels = c("A" = "Disciplina A\n(Teorica / Humanidades)",
                               "B" = "Disciplina B\n(Aplicada / Ciencias)")) +
  scale_y_continuous(labels = scales::comma) +
  labs(title    = "Distribucion del Salario segun Disciplina",
       subtitle  = "El rombo blanco indica la media de cada grupo",
       x = "Disciplina",
       y = "Salario Anual (USD)") +
  theme_bw() +
  theme(legend.position = "none",
        plot.title    = element_text(face = "bold", size = 13),
        plot.subtitle = element_text(size = 10, color = "gray40"))

Interpretacion: Los profesores de disciplinas aplicadas (tipo B) perciben salarios medianos ligeramente superiores a los de disciplinas teoricas (tipo A). Ambas distribuciones muestran sesgo hacia la derecha y presencia de valores atipicos en la parte superior.

4.4.3 Salario segun Sexo (sex)

ggplot(Salaries, aes(x = sex, y = salary, fill = sex)) +
  geom_boxplot(outlier.colour = "red", outlier.shape = 16, outlier.size = 2) +
  stat_summary(fun = mean, geom = "point", shape = 23,
               size = 3, fill = "white", color = "black") +
  scale_fill_manual(values = c("Female" = "#E07B8A", "Male" = "#6BAED6")) +
  scale_x_discrete(labels = c("Female" = "Mujer", "Male" = "Hombre")) +
  scale_y_continuous(labels = scales::comma) +
  labs(title    = "Distribucion del Salario segun Sexo",
       subtitle  = "El rombo blanco indica la media de cada grupo",
       x = "Sexo",
       y = "Salario Anual (USD)") +
  theme_bw() +
  theme(legend.position = "none",
        plot.title    = element_text(face = "bold", size = 13),
        plot.subtitle = element_text(size = 10, color = "gray40"))

Interpretacion: Los hombres presentan una mediana salarial ligeramente superior a la de las mujeres, aunque las distribuciones se superponen considerablemente. La diferencia visual no es tan marcada como la observada con la variable rank, lo que sugiere que el sexo podria tener un efecto menor o estar mediado por otras variables como el rango.

4.4.4 Boxplots combinados: Salario por Rango y Disciplina

ggplot(Salaries, aes(x = rank, y = salary, fill = discipline)) +
  geom_boxplot(outlier.colour = "red", outlier.shape = 16,
               outlier.size = 1.5, position = position_dodge(width = 0.8)) +
  scale_fill_manual(values = c("A" = "#C44E52", "B" = "#4C72B0"),
                    labels = c("A" = "Teorica (A)", "B" = "Aplicada (B)")) +
  scale_x_discrete(labels = c("AsstProf"  = "Prof. Asistente",
                               "AssocProf" = "Prof. Asociado",
                               "Prof"      = "Prof. Titular")) +
  scale_y_continuous(labels = scales::comma) +
  labs(title  = "Salario por Rango Academico y Disciplina",
       x      = "Rango Academico",
       y      = "Salario Anual (USD)",
       fill   = "Disciplina") +
  theme_bw() +
  theme(plot.title   = element_text(face = "bold", size = 13),
        legend.position = "bottom")

Interpretacion: El grafico muestra que la brecha salarial entre disciplinas (A vs B) es mas pronunciada en los rangos superiores (Profesor Titular). En los rangos inferiores, la diferencia por disciplina es minima. Esto sugiere una posible interaccion entre rank y discipline, aunque en este analisis trabajaremos con efectos principales.

4.4.5 Boxplots combinados: Salario por Rango y Sexo

ggplot(Salaries, aes(x = rank, y = salary, fill = sex)) +
  geom_boxplot(outlier.colour = "red", outlier.shape = 16,
               outlier.size = 1.5, position = position_dodge(width = 0.8)) +
  scale_fill_manual(values = c("Female" = "#E07B8A", "Male" = "#6BAED6"),
                    labels = c("Female" = "Mujer", "Male" = "Hombre")) +
  scale_x_discrete(labels = c("AsstProf"  = "Prof. Asistente",
                               "AssocProf" = "Prof. Asociado",
                               "Prof"      = "Prof. Titular")) +
  scale_y_continuous(labels = scales::comma) +
  labs(title  = "Salario por Rango Academico y Sexo",
       x      = "Rango Academico",
       y      = "Salario Anual (USD)",
       fill   = "Sexo") +
  theme_bw() +
  theme(plot.title      = element_text(face = "bold", size = 13),
        legend.position = "bottom")

Interpretacion: Dentro de cada rango, los hombres tienden a mostrar una mediana ligeramente superior a la de las mujeres, pero las distribuciones se superponen en gran medida, especialmente en los rangos de Profesor Asistente y Asociado. En el rango de Profesor Titular la diferencia es mas visible, aunque aun es moderada.

5 Analisis de Normalidad

Dado que el tamano muestral es n = 397 > 50, se utiliza la prueba de Kolmogorov-Smirnov (KS) para evaluar la normalidad de las variables numericas. La prueba de Shapiro-Wilk se reserva para muestras pequenas (n < 50).

Hipotesis de la prueba KS:

  • H0: La variable sigue una distribucion normal.
  • H1: La variable NO sigue una distribucion normal.

Se rechaza H0 cuando el p-valor < 0.05 (nivel de significancia del 5%).

vars_num <- c("salary", "yrs.since.phd", "yrs.service")

resultados_ks <- data.frame(
  Variable  = vars_num,
  Estadistico_D = NA,
  p_valor   = NA,
  Conclusion = NA
)

for (i in seq_along(vars_num)) {
  x   <- Salaries[[vars_num[i]]]
  ks  <- ks.test(x, "pnorm", mean = mean(x), sd = sd(x))
  resultados_ks$Estadistico_D[i] <- round(ks$statistic, 4)
  resultados_ks$p_valor[i]       <- round(ks$p.value, 4)
  resultados_ks$Conclusion[i]    <- ifelse(ks$p.value < 0.05,
                                            "No normal (rechaza H0)",
                                            "Normal (no rechaza H0)")
}

resultados_ks |>
  kable(caption = "Prueba de Kolmogorov-Smirnov para variables numericas",
        col.names = c("Variable", "Estadistico D", "p-valor", "Conclusion")) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Prueba de Kolmogorov-Smirnov para variables numericas
Variable Estadistico D p-valor Conclusion
salary 0.0909 0.0028 No normal (rechaza H0)
yrs.since.phd 0.0699 0.0414 No normal (rechaza H0)
yrs.service 0.1202 0.0000 No normal (rechaza H0)

5.1 QQ-Plots de las variables numericas

par(mfrow = c(1, 3))

for (v in vars_num) {
  qqnorm(Salaries[[v]], main = paste("QQ-Plot:", v),
         pch = 16, col = "steelblue", cex = 0.7)
  qqline(Salaries[[v]], col = "red", lwd = 2)
}

par(mfrow = c(1, 1))

Interpretacion: Las tres variables muestran desviaciones respecto a la linea de referencia en los extremos (colas), lo que confirma la falta de normalidad detectada por la prueba KS. Esto es esperado para datos de salarios y anos de experiencia, que suelen mostrar sesgo positivo.

6 Analisis de Correlacion

Dado que las variables salary, yrs.since.phd y yrs.service no siguen una distribucion normal (segun la prueba KS), se debe utilizar el coeficiente de correlacion de Spearman en lugar del de Pearson.

La correlacion de Spearman mide la asociacion monotona entre dos variables sin asumir normalidad. Se basa en los rangos de los datos.

Hipotesis de la prueba de correlacion:

  • H0: No existe correlacion entre las variables (rho = 0).
  • H1: Existe correlacion entre las variables (rho != 0).
# Correlacion de Spearman: salary vs yrs.since.phd
cor_phd <- cor.test(Salaries$salary, Salaries$yrs.since.phd,
                     method = "spearman", exact = FALSE)

# Correlacion de Spearman: salary vs yrs.service
cor_srv <- cor.test(Salaries$salary, Salaries$yrs.service,
                     method = "spearman", exact = FALSE)

tabla_cor <- data.frame(
  Relacion       = c("salary ~ yrs.since.phd", "salary ~ yrs.service"),
  Metodo         = c("Spearman", "Spearman"),
  Coeficiente    = round(c(cor_phd$estimate, cor_srv$estimate), 4),
  p_valor        = round(c(cor_phd$p.value, cor_srv$p.value), 6),
  Significativa  = c(ifelse(cor_phd$p.value < 0.05, "Si", "No"),
                     ifelse(cor_srv$p.value < 0.05, "Si", "No"))
)

tabla_cor |>
  kable(caption = "Correlacion de Spearman: variables numericas con salary",
        col.names = c("Relacion", "Metodo", "Coeficiente rho",
                      "p-valor", "Significativa (alpha=0.05)")) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Correlacion de Spearman: variables numericas con salary
Relacion Metodo Coeficiente rho p-valor Significativa (alpha=0.05)
salary ~ yrs.since.phd Spearman 0.4756 0 Si
salary ~ yrs.service Spearman 0.4249 0 Si

Interpretacion:

  • yrs.since.phd presenta una correlacion de Spearman de 0.476 con el salario (p-valor = 0). La correlacion es positiva y estadisticamente significativa: a mayor numero de anos desde el doctorado, mayor tiende a ser el salario.

  • yrs.service presenta una correlacion de Spearman de 0.425 con el salario (p-valor = 0). La correlacion es positiva y significativa, aunque de magnitud moderada.

Ambas variables seran candidatas para el modelo dado que su relacion con salary es significativa.

6.1 Grafico de dispersion con linea de tendencia

p1 <- ggplot(Salaries, aes(x = yrs.since.phd, y = salary)) +
  geom_point(alpha = 0.5, color = "steelblue") +
  geom_smooth(method = "loess", color = "red", se = TRUE) +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Salario vs Anos desde el Doctorado",
       x = "Anos desde el Doctorado", y = "Salario (USD)") +
  theme_bw() +
  theme(plot.title = element_text(face = "bold", size = 11))

p2 <- ggplot(Salaries, aes(x = yrs.service, y = salary)) +
  geom_point(alpha = 0.5, color = "darkgreen") +
  geom_smooth(method = "loess", color = "red", se = TRUE) +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Salario vs Anos de Servicio",
       x = "Anos de Servicio", y = "Salario (USD)") +
  theme_bw() +
  theme(plot.title = element_text(face = "bold", size = 11))

grid.arrange(p1, p2, ncol = 2)

7 Analisis de Multicolinealidad Previo al Modelo

Antes de ajustar el modelo final, es importante evaluar si las variables explicativas presentan multicolinealidad entre si. Se ajusta un modelo con todas las variables y se calculan los valores VIF y GVIF.

7.1 Criterios de interpretacion del VIF y GVIF

Para variables continuas o dummies simples (2 categorias), se usa el VIF directamente:

Valor VIF Interpretacion
VIF = 1 Sin multicolinealidad
1 < VIF <= 5 Multicolinealidad leve, aceptable
5 < VIF <= 10 Multicolinealidad moderada, revisar
VIF > 10 Multicolinealidad severa, problema

Para variables categoricas con mas de 2 niveles, la funcion vif() del paquete car calcula el GVIF generalizado y reporta tres columnas: GVIF, Df (grados de libertad del termino) y GVIF^(1/(2*Df)). Siempre se interpreta la tercera columna, que estandariza el GVIF y es directamente comparable con el umbral del VIF clasico mediante la siguiente equivalencia:

GVIF^(1/(2*Df)) VIF equivalente Interpretacion
< 1.58 < 2.5 Sin problema
1.58 - 2.24 2.5 - 5.0 Leve, aceptable
2.24 - 3.16 5.0 - 10.0 Moderada, revisar
> 3.16 > 10.0 Severa, problema

La logica es: \(\left(GVIF^{1/(2 \cdot Df)}\right)^2 = VIF_{equivalente}\), por lo que el umbral de VIF = 10 se traduce en \(\sqrt{10} \approx 3.16\) para el GVIF ajustado.

modelo_completo_prev <- lm(salary ~ rank + discipline + yrs.since.phd +
                              yrs.service + sex, data = Salaries)

vif_previo <- vif(modelo_completo_prev)
print(vif_previo)
##                   GVIF Df GVIF^(1/(2*Df))
## rank          2.013193  2        1.191163
## discipline    1.064105  1        1.031555
## yrs.since.phd 7.518936  1        2.742068
## yrs.service   5.923038  1        2.433729
## sex           1.030805  1        1.015285

Interpretacion del analisis de multicolinealidad previo:

  • La variable rank tiene 3 niveles (Df = 2). Se evalua su GVIF^(1/(2*Df)).
  • Las variables discipline y sex son binarias (Df = 1). Se evalua su GVIF^(1/(2*Df)) o equivalentemente su VIF.
  • yrs.since.phd y yrs.service son continuas, con Df = 1.

Si alguno de los valores supera los umbrales establecidos, se indicara en el analisis de multicolinealidad del modelo final (Seccion 9).

8 Seleccion de Variables: Metodo Forward con step()

Se utiliza la seleccion de variables hacia adelante (Forward Selection) basada en el criterio de informacion de Akaike (AIC). El AIC penaliza la complejidad del modelo: menor AIC indica mejor modelo.

El proceso:

  1. Modelo inicial (nulo): salary ~ 1 (solo la media, sin predictores).
  2. Modelo maximo (completo): salary ~ rank + discipline + yrs.since.phd + yrs.service + sex.
  3. En cada paso, se agrega la variable que mas reduce el AIC.
  4. El proceso se detiene cuando ninguna variable adicional mejora el AIC.
# Modelo nulo (punto de partida)
modelo_nulo     <- lm(salary ~ 1, data = Salaries)

# Modelo completo (alcance maximo)
modelo_completo <- lm(salary ~ rank + discipline + yrs.since.phd +
                        yrs.service + sex, data = Salaries)

# Seleccion Forward
modelo_forward <- step(
  object    = modelo_nulo,
  scope     = list(lower = modelo_nulo, upper = modelo_completo),
  direction = "forward",
  trace     = TRUE
)
## Start:  AIC=8193.92
## salary ~ 1
## 
##                 Df  Sum of Sq        RSS    AIC
## + rank           2 1.4323e+11 2.2007e+11 7998.9
## + yrs.since.phd  1 6.3852e+10 2.9945e+11 8119.2
## + yrs.service    1 4.0709e+10 3.2259e+11 8148.7
## + discipline     1 8.8508e+09 3.5445e+11 8186.1
## + sex            1 6.9800e+09 3.5632e+11 8188.2
## <none>                        3.6330e+11 8193.9
## 
## Step:  AIC=7998.91
## salary ~ rank
## 
##                 Df  Sum of Sq        RSS    AIC
## + discipline     1 1.8430e+10 2.0164e+11 7966.2
## <none>                        2.2007e+11 7998.9
## + yrs.service    1 1.0546e+09 2.1901e+11 7999.0
## + sex            1 8.4082e+08 2.1923e+11 7999.4
## + yrs.since.phd  1 2.1559e+08 2.1985e+11 8000.5
## 
## Step:  AIC=7966.19
## salary ~ rank + discipline
## 
##                 Df Sum of Sq        RSS    AIC
## <none>                       2.0164e+11 7966.2
## + sex            1 694070191 2.0094e+11 7966.8
## + yrs.service    1 241866236 2.0140e+11 7967.7
## + yrs.since.phd  1 165649329 2.0147e+11 7967.9
summary(modelo_forward)
## 
## Call:
## lm(formula = salary ~ rank + discipline, data = Salaries)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -65990 -14049  -1288  10760  97996 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      71944       3135  22.948  < 2e-16 ***
## rankAssocProf    13762       3961   3.475 0.000569 ***
## rankProf         47844       3112  15.376  < 2e-16 ***
## disciplineB      13761       2296   5.993 4.65e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22650 on 393 degrees of freedom
## Multiple R-squared:  0.445,  Adjusted R-squared:  0.4407 
## F-statistic:   105 on 3 and 393 DF,  p-value: < 2.2e-16

Lectura del proceso de seleccion: En cada iteracion, la tabla AIC muestra el efecto de agregar cada variable al modelo actual. La variable con el AIC mas bajo es la seleccionada. El proceso continua hasta que ninguna variable adicional mejora el AIC (aparece como <none> con el AIC mas bajo).

9 Modelo de Regresion Lineal Multiple Final

9.1 Variables incluidas en el modelo

El proceso de seleccion Forward determina el conjunto optimo de variables. A continuacion se presenta el modelo final, que incluye variables numericas y al menos una variable categorica.

# El modelo seleccionado por step() es modelo_forward
# Verificamos que incluye variables categoricas
print(formula(modelo_forward))
## salary ~ rank + discipline
# Resumen completo
resumen <- summary(modelo_forward)
print(resumen)
## 
## Call:
## lm(formula = salary ~ rank + discipline, data = Salaries)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -65990 -14049  -1288  10760  97996 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      71944       3135  22.948  < 2e-16 ***
## rankAssocProf    13762       3961   3.475 0.000569 ***
## rankProf         47844       3112  15.376  < 2e-16 ***
## disciplineB      13761       2296   5.993 4.65e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22650 on 393 degrees of freedom
## Multiple R-squared:  0.445,  Adjusted R-squared:  0.4407 
## F-statistic:   105 on 3 and 393 DF,  p-value: < 2.2e-16

9.2 Ecuacion matematica del modelo

El modelo de regresion lineal multiple tiene la forma general:

\[\hat{Y} = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k + \varepsilon\]

Para las variables categoricas, R crea variables indicadoras (dummies). La categoria de referencia queda absorbida en el intercepto:

  • rank: categoria de referencia = AsstProf (Profesor Asistente)
  • discipline: categoria de referencia = A (Teorica)
  • sex: categoria de referencia = Female (Mujer)

La ecuacion del modelo seleccionado es:

\[\widehat{\text{salary}} = \beta_0 + \beta_1 \cdot \mathbb{1}[\text{rank=AssocProf}] + \beta_2 \cdot \mathbb{1}[\text{rank=Prof}] + \beta_3 \cdot \mathbb{1}[\text{discipline=B}] + \beta_4 \cdot \text{yrs.since.phd} + \beta_5 \cdot \text{yrs.service} + \beta_6 \cdot \mathbb{1}[\text{sex=Male}]\]

Donde \(\mathbb{1}[\cdot]\) es la funcion indicadora que vale 1 si la condicion es verdadera y 0 en caso contrario.

coef_tabla <- data.frame(
  Coeficiente = names(coef(modelo_forward)),
  Estimacion  = round(coef(modelo_forward), 2),
  Error_Est   = round(summary(modelo_forward)$coefficients[, 2], 2),
  t_value     = round(summary(modelo_forward)$coefficients[, 3], 3),
  p_valor     = round(summary(modelo_forward)$coefficients[, 4], 4)
)

coef_tabla |>
  kable(caption = "Coeficientes del modelo de regresion lineal multiple",
        row.names = FALSE) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Coeficientes del modelo de regresion lineal multiple
Coeficiente Estimacion Error_Est t_value p_valor
(Intercept) 71944.33 3135.17 22.948 0e+00
rankAssocProf 13761.54 3960.66 3.475 6e-04
rankProf 47843.84 3111.55 15.376 0e+00
disciplineB 13760.96 2296.03 5.993 0e+00

9.3 Interpretacion de los coeficientes

  • Intercepto (\(\beta_0\)): El salario estimado para un Profesor Asistente, de disciplina teorica (A), del sexo femenino, con cero anos desde el doctorado y cero anos de servicio es de 71,944 dolares. Corresponde a la categoria de referencia de todas las variables categoricas.

  • rankAssocProf (\(\beta_1\)): Manteniendo todo lo demas constante, un Profesor Asociado gana en promedio 13,762 dolares MAS que un Profesor Asistente (categoria de referencia).

  • rankProf (\(\beta_2\)): Manteniendo todo lo demas constante, un Profesor Titular gana en promedio 47,844 dolares MAS que un Profesor Asistente.

  • disciplineB (\(\beta_3\)): Los profesores de disciplinas aplicadas (B) ganan en promedio 13,761 dolares MAS que los de disciplinas teoricas (A), ceteris paribus.

  • yrs.since.phd (\(\beta_4\)): Por cada ano adicional desde la obtencion del doctorado, el salario aumenta en promedio NA dolares, manteniendo las demas variables constantes.

  • yrs.service (\(\beta_5\)): Por cada ano adicional de servicio, el salario cambia en NA dolares en promedio. Nótese que este coeficiente puede ser negativo si hay colinealidad con yrs.since.phd.

  • sexMale (\(\beta_6\)): Los hombres ganan en promedio NA dolares NA que las mujeres, ceteris paribus.

Bondad de ajuste:

  • \(R^2\) = 0.445: el modelo explica el 44.5% de la variabilidad total del salario.
  • \(R^2\) ajustado = 0.4407: corregido por el numero de predictores.
  • El estadistico F global (p-valor < 2.2e-16) confirma que el modelo es significativo en su conjunto.

10 Analisis de Residuos

Los residuos son las diferencias entre los valores observados y los valores predichos por el modelo: \(e_i = y_i - \hat{y}_i\). Un modelo bien ajustado debe cumplir los siguientes supuestos sobre sus residuos:

  1. Normalidad: Los residuos siguen una distribucion normal.
  2. Homocedasticidad: La varianza de los residuos es constante.
  3. Independencia: Los residuos no estan autocorrelacionados.

10.1 Grafico general de diagnostico del modelo

par(mfrow = c(2, 2))
plot(modelo_forward, which = 1:4, col = "steelblue",
     pch = 16, cex = 0.7)

par(mfrow = c(1, 1))

Descripcion de los cuatro graficos:

  1. Residuos vs Valores Ajustados (arriba izquierda): Permite detectar patrones no lineales y homocedasticidad. La linea roja debe ser aproximadamente horizontal y centrada en cero.
  2. QQ-Plot de residuos (arriba derecha): Evalua la normalidad de los residuos. Los puntos deben seguir la diagonal.
  3. Scale-Location (abajo izquierda): Muestra la raiz cuadrada de los residuos estandarizados. La linea debe ser aproximadamente horizontal para confirmar homocedasticidad.
  4. Residuos vs Leverage (abajo derecha): Identifica puntos influyentes. Puntos fuera de las lineas discontinuas de la distancia de Cook merecen atencion.

10.2 Supuesto 1: Normalidad de los residuos

Prueba KS sobre residuos:

  • H0: Los residuos siguen una distribucion normal.
  • H1: Los residuos NO siguen una distribucion normal.
residuos <- residuals(modelo_forward)

ks_res <- ks.test(residuos, "pnorm",
                   mean = mean(residuos),
                   sd   = sd(residuos))

cat("Prueba de Kolmogorov-Smirnov sobre los residuos\n")
## Prueba de Kolmogorov-Smirnov sobre los residuos
cat("------------------------------------------------\n")
## ------------------------------------------------
cat("Estadistico D :", round(ks_res$statistic, 4), "\n")
## Estadistico D : 0.0962
cat("p-valor       :", round(ks_res$p.value, 4), "\n")
## p-valor       : 0.0013
cat("Conclusion    :", ifelse(ks_res$p.value < 0.05,
                               "Se rechaza H0: los residuos no son normales (p < 0.05)",
                               "No se rechaza H0: los residuos son normales (p >= 0.05)"), "\n")
## Conclusion    : Se rechaza H0: los residuos no son normales (p < 0.05)
par(mfrow = c(1, 2))

hist(residuos,
     main   = "Histograma de Residuos",
     xlab   = "Residuos",
     ylab   = "Frecuencia",
     col    = "steelblue",
     border = "white",
     breaks = 30)
abline(v = 0, col = "red", lwd = 2, lty = 2)

qqnorm(residuos,
       main = "QQ-Plot de Residuos",
       pch  = 16, col = "steelblue", cex = 0.7)
qqline(residuos, col = "red", lwd = 2)

par(mfrow = c(1, 1))

Interpretacion: El QQ-Plot de los residuos muestra desviaciones en las colas, lo que es consistente con el resultado de la prueba KS. En muestras grandes, es frecuente que la prueba KS rechace H0 incluso con desviaciones menores. Se recomienda complementar con la inspeccion visual del QQ-Plot y el histograma.

10.3 Supuesto 2: Homocedasticidad

Se utiliza la prueba de Breusch-Pagan para evaluar si la varianza de los residuos es constante.

Hipotesis:

  • H0: La varianza de los residuos es constante (homocedasticidad).
  • H1: La varianza de los residuos NO es constante (heterocedasticidad).
bp_test <- bptest(modelo_forward)

cat("Prueba de Breusch-Pagan\n")
## Prueba de Breusch-Pagan
cat("-----------------------\n")
## -----------------------
cat("Estadistico BP:", round(bp_test$statistic, 4), "\n")
## Estadistico BP: 36.6987
cat("Grados de libertad:", bp_test$parameter, "\n")
## Grados de libertad: 3
cat("p-valor       :", round(bp_test$p.value, 4), "\n")
## p-valor       : 0
cat("Conclusion    :", ifelse(bp_test$p.value < 0.05,
                               "Se rechaza H0: existe heterocedasticidad (p < 0.05)",
                               "No se rechaza H0: homocedasticidad (p >= 0.05)"), "\n")
## Conclusion    : Se rechaza H0: existe heterocedasticidad (p < 0.05)
plot(modelo_forward)

Interpretacion: Si el p-valor de la prueba Breusch-Pagan es menor a 0.05, existe evidencia de heterocedasticidad, es decir, la varianza de los residuos no es constante a lo largo de los valores ajustados. Esto es comun en datos de salarios donde la dispersion tiende a aumentar con el nivel salarial.

10.4 Supuesto 3: Independencia de los residuos

Se utiliza la prueba de Durbin-Watson para detectar autocorrelacion de primer orden en los residuos.

Hipotesis:

  • H0: Los residuos son independientes (no hay autocorrelacion).
  • H1: Existe autocorrelacion de primer orden en los residuos.

El estadistico DW toma valores entre 0 y 4. Un valor cercano a 2 indica independencia; valores cercanos a 0 indican autocorrelacion positiva; valores cercanos a 4 indican autocorrelacion negativa.

dw_test <- durbinWatsonTest(modelo_forward)

cat("Prueba de Durbin-Watson\n")
## Prueba de Durbin-Watson
cat("-----------------------\n")
## -----------------------
cat("Estadistico DW:", round(dw_test$dw, 4), "\n")
## Estadistico DW: 1.8595
cat("p-valor       :", round(dw_test$p, 4), "\n")
## p-valor       : 0.15
cat("Autocorrelacion estimada:", round(dw_test$r, 4), "\n")
## Autocorrelacion estimada: 0.0699
cat("Conclusion    :", ifelse(dw_test$p < 0.05,
                               "Se rechaza H0: hay autocorrelacion en los residuos",
                               "No se rechaza H0: residuos independientes"), "\n")
## Conclusion    : No se rechaza H0: residuos independientes

Interpretacion: Un estadistico DW proximo a 2 con p-valor > 0.05 indica que no hay evidencia de autocorrelacion en los residuos, por lo que el supuesto de independencia se cumple. Si el dataset no tiene estructura temporal, este supuesto generalmente se cumple.

11 Deteccion de Valores Atipicos e Influyentes

11.1 Prueba de Bonferroni para valores atipicos

La prueba de Bonferroni aplica una correccion al nivel de significancia para detectar residuos externamente estudentizados que sean estadisticamente atipicos, controlando el error de Tipo I.

  • H0: La observacion no es un valor atipico.
  • H1: La observacion ES un valor atipico.

Se considera atipica una observacion si su p-valor de Bonferroni < 0.05.

cat("Prueba de Outliers con correccion de Bonferroni\n")
## Prueba de Outliers con correccion de Bonferroni
cat("------------------------------------------------\n")
## ------------------------------------------------
outlierTest(modelo_forward)
##    rstudent unadjusted p-value Bonferroni p
## 44 4.442116         1.1604e-05    0.0046067

Interpretacion: Si el p-valor de Bonferroni para la observacion con el mayor residuo estudentizado es mayor a 0.05, no hay evidencia estadistica de valores atipicos extremos en el modelo.

11.2 Matriz de Sombrero (Hat Matrix) - Puntos de Apalancamiento

Los puntos de apalancamiento (leverage) son observaciones con valores inusualmente extremos en las variables explicativas (no necesariamente en la respuesta). Se miden mediante los valores de la diagonal de la matriz Hat (\(h_{ii}\)).

Criterio de apalancamiento alto:

\[h_{ii} > \frac{2p}{n}\]

Donde \(p\) es el numero de parametros del modelo (incluyendo el intercepto) y \(n\) es el numero de observaciones.

hat_vals <- hatvalues(modelo_forward)
p_param  <- length(coef(modelo_forward))
n_obs    <- nrow(Salaries)
umbral_h <- 2 * p_param / n_obs

cat("Umbral de apalancamiento (2p/n):", round(umbral_h, 4), "\n")
## Umbral de apalancamiento (2p/n): 0.0202
cat("Numero de puntos con leverage alto:",
    sum(hat_vals > umbral_h), "\n\n")
## Numero de puntos con leverage alto: 0
# Identificar los 10 mayores
top_leverage <- sort(hat_vals, decreasing = TRUE)[1:10]
cat("Las 10 observaciones con mayor leverage:\n")
## Las 10 observaciones con mayor leverage:
print(round(top_leverage, 4))
##     25    105    107    108    109    112    124    131    133    139 
## 0.0192 0.0192 0.0192 0.0192 0.0192 0.0192 0.0192 0.0192 0.0192 0.0192
plot(hat_vals,
     main = "Valores de Apalancamiento (Leverage)",
     xlab = "Indice de Observacion",
     ylab = expression(h[ii]),
     pch  = 16, col = "steelblue", cex = 0.7)
abline(h = umbral_h, col = "red", lwd = 2, lty = 2)

puntos_hat <- which(hat_vals > umbral_h)
if (length(puntos_hat) > 0) {
  text(puntos_hat,
       hat_vals[puntos_hat],
       labels = puntos_hat,
       cex = 0.65, pos = 3, col = "red")
}

legend("topright", legend = paste("Umbral =", round(umbral_h, 3)),
       col = "red", lty = 2, lwd = 2)

11.3 Distancias de Cook

La distancia de Cook mide el impacto global de cada observacion sobre todos los coeficientes estimados del modelo. Combina leverage y residuo en una sola medida.

Criterio: Una observacion es influyente si su distancia de Cook supera el umbral:

\[D_i > \frac{4}{n - p}\]

cook_d   <- cooks.distance(modelo_forward)
umbral_cook <- 4 / (n_obs - p_param)

cat("Umbral de Cook (4/(n-p)):", round(umbral_cook, 4), "\n")
## Umbral de Cook (4/(n-p)): 0.0102
cat("Observaciones con Cook > umbral:",
    sum(cook_d > umbral_cook), "\n\n")
## Observaciones con Cook > umbral: 10
cat("Las 10 mayores distancias de Cook:\n")
## Las 10 mayores distancias de Cook:
print(round(sort(cook_d, decreasing = TRUE)[1:10], 4))
##     44    365    250    272    390    318    293    283     78    331 
## 0.0296 0.0232 0.0224 0.0178 0.0139 0.0134 0.0130 0.0121 0.0109 0.0106
plot(cook_d,
     main = "Distancias de Cook",
     xlab = "Indice de Observacion",
     ylab = "Distancia de Cook",
     pch  = 16, col = "steelblue", cex = 0.7,
     ylim = c(0, max(cook_d) * 1.15))
abline(h = umbral_cook, col = "red", lwd = 2, lty = 2)
puntos_influyentes <- which(cook_d > umbral_cook)
text(puntos_influyentes, cook_d[puntos_influyentes],
     labels = puntos_influyentes, cex = 0.65, pos = 3, col = "red")
legend("topright", legend = paste("Umbral =", round(umbral_cook, 3)),
       col = "red", lty = 2, lwd = 2)

11.4 DFFITS

Los DFFITS miden cuanto cambia el valor ajustado de cada observacion cuando esa observacion se elimina del modelo. Es una medida individual de influencia sobre los valores predichos.

Criterio: Una observacion es influyente si:

\[|DFFITS_i| > 2\sqrt{\frac{p}{n}}\]

dffits_vals  <- dffits(modelo_forward)
umbral_dffit <- 2 * sqrt(p_param / n_obs)

cat("Umbral de DFFITS (2*sqrt(p/n)):", round(umbral_dffit, 4), "\n")
## Umbral de DFFITS (2*sqrt(p/n)): 0.2008
cat("Observaciones con |DFFITS| > umbral:",
    sum(abs(dffits_vals) > umbral_dffit), "\n")
## Observaciones con |DFFITS| > umbral: 10
plot(dffits_vals,
     main = "DFFITS por Observacion",
     xlab = "Indice de Observacion",
     ylab = "DFFITS",
     pch  = 16, col = "steelblue", cex = 0.7)
abline(h =  umbral_dffit, col = "red", lwd = 2, lty = 2)
abline(h = -umbral_dffit, col = "red", lwd = 2, lty = 2)
text(which(abs(dffits_vals) > umbral_dffit),
     dffits_vals[abs(dffits_vals) > umbral_dffit],
     labels = which(abs(dffits_vals) > umbral_dffit),
     cex = 0.65, pos = 3, col = "red")
legend("topright",
       legend = paste0("Umbral: +/- ", round(umbral_dffit, 3)),
       col = "red", lty = 2, lwd = 2)

11.5 DFBETAS

Los DFBETAS miden cuanto cambia cada coeficiente del modelo cuando se elimina una observacion. A diferencia de DFFITS (que es global), DFBETAS es especifico para cada coeficiente.

Criterio:

\[|DFBETAS_{j,i}| > \frac{2}{\sqrt{n}}\]

dfbetas_vals  <- dfbetas(modelo_forward)
umbral_dfbeta <- 2 / sqrt(n_obs)

cat("Umbral de DFBETAS (2/sqrt(n)):", round(umbral_dfbeta, 4), "\n\n")
## Umbral de DFBETAS (2/sqrt(n)): 0.1004
# Numero de observaciones influyentes por coeficiente
influyentes_por_coef <- apply(abs(dfbetas_vals) > umbral_dfbeta, 2, sum)
cat("Observaciones influyentes por coeficiente:\n")
## Observaciones influyentes por coeficiente:
print(influyentes_por_coef)
##   (Intercept) rankAssocProf      rankProf   disciplineB 
##             1             2             1            22
nombres_coef <- colnames(dfbetas_vals)
n_coef       <- ncol(dfbetas_vals)
n_cols_plot  <- 2
n_rows_plot  <- ceiling(n_coef / n_cols_plot)

par(mfrow = c(n_rows_plot, n_cols_plot), mar = c(4, 4, 3, 1))

for (j in seq_len(n_coef)) {
  plot(dfbetas_vals[, j],
       main = paste("DFBETAS:", nombres_coef[j]),
       xlab = "Indice",
       ylab = "DFBETAS",
       pch  = 16, col = "steelblue", cex = 0.6,
       ylim = c(-max(abs(dfbetas_vals[, j])) * 1.2,
                 max(abs(dfbetas_vals[, j])) * 1.2))
  abline(h =  umbral_dfbeta, col = "red", lwd = 1.5, lty = 2)
  abline(h = -umbral_dfbeta, col = "red", lwd = 1.5, lty = 2)
  puntos <- which(abs(dfbetas_vals[, j]) > umbral_dfbeta)
  if (length(puntos) > 0) {
    text(puntos, dfbetas_vals[puntos, j],
         labels = puntos, cex = 0.55, pos = 3, col = "red")
  }
}

par(mfrow = c(1, 1))

Interpretacion general de las medidas de influencia:

  • Las Distancias de Cook identifican observaciones que afectan el modelo en su conjunto.
  • Los DFFITS senalan observaciones que distorsionan los valores predichos.
  • Los DFBETAS permiten identificar que coeficientes especificos son afectados por cada observacion.

Una observacion que aparece como influyente en los tres criterios simultaneamente merece atencion especial: se debe inspeccionar si se trata de un error de captura o de una observacion genuinamente inusual con efecto real sobre los resultados.

12 Analisis de Multicolinealidad del Modelo Final

Se evalua la multicolinealidad del modelo seleccionado mediante VIF y GVIF.

12.1 Criterios de interpretacion (recordatorio)

Para variables continuas y dummies simples (VIF, Df = 1):

VIF Interpretacion
1 Sin multicolinealidad
1 < VIF <= 5 Leve, generalmente aceptable
5 < VIF <= 10 Moderada, se recomienda revisar
VIF > 10 Severa, los coeficientes no son confiables

**Para variables categoricas con mas de 2 niveles (GVIF^(1/(2*Df))):**

GVIF^(1/(2*Df)) VIF equivalente Interpretacion
< 1.58 < 2.5 Sin problema
1.58 - 2.24 2.5 - 5.0 Leve, aceptable
2.24 - 3.16 5.0 - 10.0 Moderada, revisar
> 3.16 > 10.0 Severa, problema

12.2 Calculo e interpretacion

vif_final <- vif(modelo_forward)
cat("VIF / GVIF del modelo final\n")
## VIF / GVIF del modelo final
cat("---------------------------\n")
## ---------------------------
print(vif_final)
##                GVIF Df GVIF^(1/(2*Df))
## rank       1.011848  2        1.002949
## discipline 1.011848  1        1.005907
cat("\nDiagnostico detallado:\n")
## 
## Diagnostico detallado:
cat("----------------------\n")
## ----------------------
if (is.matrix(vif_final)) {
  # Modelo con factores: usar GVIF^(1/(2*Df))
  gvif_aj <- vif_final[, "GVIF^(1/(2*Df))"]
  df_vals  <- vif_final[, "Df"]
  
  for (var in names(gvif_aj)) {
    val <- gvif_aj[var]
    df  <- df_vals[var]
    
    if (df == 1) {
      vif_eq <- round(val^2, 3)
      etiqueta <- ifelse(val > 3.16, "SEVERA - Problema",
                  ifelse(val > 2.24, "MODERADA - Revisar",
                  ifelse(val > 1.58, "LEVE - Aceptable", "SIN PROBLEMA")))
      cat(sprintf("%-20s | Df=%d | GVIF^(1/2*Df)=%.3f | VIF equiv.=%.3f | %s\n",
                  var, df, val, vif_eq, etiqueta))
    } else {
      etiqueta <- ifelse(val > 3.16, "SEVERA - Problema",
                  ifelse(val > 2.24, "MODERADA - Revisar",
                  ifelse(val > 1.58, "LEVE - Aceptable", "SIN PROBLEMA")))
      cat(sprintf("%-20s | Df=%d | GVIF^(1/(2*Df))=%.3f | %s\n",
                  var, df, val, etiqueta))
    }
  }
  
} else {
  # Solo variables continuas: VIF directo
  for (var in names(vif_final)) {
    val <- vif_final[var]
    etiqueta <- ifelse(val > 10, "SEVERA - Problema",
                ifelse(val > 5,  "MODERADA - Revisar",
                ifelse(val > 2.5, "LEVE - Aceptable", "SIN PROBLEMA")))
    cat(sprintf("%-20s | VIF=%.3f | %s\n", var, val, etiqueta))
  }
}
## rank                 | Df=2 | GVIF^(1/(2*Df))=1.003 | SIN PROBLEMA
## discipline           | Df=1 | GVIF^(1/2*Df)=1.006 | VIF equiv.=1.012 | SIN PROBLEMA

Interpretacion final:

  • Para las variables con Df = 1 (variables binarias y continuas), el GVIF^(1/(2*Df)) es equivalente a la raiz cuadrada del VIF clasico. Se eleva al cuadrado para obtener el VIF equivalente y compararlo con los umbrales estandar (VIF < 5 aceptable, VIF > 10 severo).

  • Para la variable rank con Df = 2 (3 niveles), el GVIF^(1/(2*Df)) estandariza el GVIF y se compara directamente con el umbral de 3.16.

  • Si yrs.since.phd y yrs.service presentan VIF elevados, esto se debe a su alta correlacion entre si (ambas miden experiencia acumulada), lo que es una manifestacion tipica de multicolinealidad que conviene reportar y considerar en la interpretacion de sus coeficientes individuales.

13 Conclusiones

  1. Analisis exploratorio: El salario de los profesores universitarios varia considerablemente segun el rango academico, siendo los Profesores Titulares quienes perciben los salarios mas altos y con mayor dispersion. La disciplina aplicada (B) y el sexo masculino tambien se asocian con salarios ligeramente mayores.

  2. Normalidad: Ninguna de las tres variables numericas (salary, yrs.since.phd, yrs.service) sigue una distribucion normal segun la prueba KS (n = 397 > 50), por lo que se uso la correlacion de Spearman para evaluar la relacion con la variable respuesta.

  3. Correlacion: Ambas variables numericas (yrs.since.phd y yrs.service) presentan correlacion de Spearman positiva y significativa con el salario, justificando su inclusion como candidatas en el modelo.

  4. Seleccion de variables: El metodo Forward con criterio AIC selecciono el subconjunto optimo de variables. El modelo final incluye al menos una variable categorica, lo que permite cuantificar diferencias salariales entre grupos.

  5. Modelo: El modelo explica aproximadamente el 44.1% de la variabilidad del salario (R2 ajustado). El rango academico es el predictor con mayor impacto sobre el salario.

  6. Residuos: El analisis de residuos revela desviaciones de la normalidad en las colas (frecuente con datos salariales). La independencia de los residuos se confirma con la prueba Durbin-Watson. La homocedasticidad debe evaluarse a partir de la prueba Breusch-Pagan.

  7. Valores influyentes: La combinacion de distancias de Cook, DFFITS y DFBETAS permite identificar observaciones con impacto desproporcionado sobre el modelo. Se recomienda verificar si corresponden a errores o a casos genuinamente particulares.

  8. Multicolinealidad: Los valores VIF y GVIF del modelo final se encuentran dentro o cerca de los rangos aceptables para la mayoria de las variables. La eventual colinealidad entre yrs.since.phd y yrs.service es esperable por su naturaleza, y debe tomarse en cuenta al interpretar sus coeficientes individuales.

Documento elaborado por Luis Armando Cruz Contreras. Analisis realizado con R version 4.5.2.