Voy a enfocar la respuesta en ¿Por qué R da un resultado "valido" cuando la teoría estadística indica que la aproximación es incorrecta?
Creo que es relevante tomando como punto de partida que R es un lenguaje de programación de dominio específico y que todavía es apropiado entenderlo como un lenguaje hecho por estadísticos para estadísticos, aunque ya haga muchas cosas más.
¿Por qué obtienes un output en lugar de un mensaje de error?
Porque la función lm()
solo evalúa si lo que defines como y en tu fórmula pertenece a la clase numeric
. Si es así procede al ajuste del modelo, ya que matemáticamente es posible resolverlo. lm()
no tiene forma de saber de dónde salieron esos números y si están cumpliendo o no los supuestos de un modelo lineal.
En cambio si y es de la clase factor y usas lm()
obtendrás un error. ¿Por qué? Porque en ese caso lm()
tiene información para identificar que se está violando uno de los supuestos de los modelos lineales. Quienes crearon la función lm()
(Chambers en S adaptada a R por Ihaka) son estadísticos y entonces les preocupa más la consistencia en la estimación de los errores estándar, que es muy vulnerable a las violaciones de supuestos, que la capacidad predictiva de un modelo.
¿Es válido el resultado en términos estadísticos?
En tu caso definitivamente no y no hay manera de que lo sea. El problema es que al asignar números a tus categorías (tal que A=1 ... D=5) introduces, de pasada, una magnitud: 1 es la mitad de 2, 4 el doble que 2, etc. El modelo lineal que ajusta lm()
va a usar esa información, porque lo que estima es E(y|x) y asume que la expectativa de y es la media de y. Y cuando se calcula una media se consideran las magnitudes (sum(x)/n). Ahora bien, si no puedes afirmar que A es la mitad que B y que E es el doble que B entonces los resultados no tienen sentido matemático ni sustantivo.
Hay un caso especial que es diferente al que planteas y es el de y dicotómica. Si y admite solo dos categorías (digamos: Sí o No) entonces puedes codificarlas como Sí = 0
y No = 1
. La clave es que la media aritmética de una variable así codificada es la proporción de y=1, ya que cuando haces la sumatoria para la media estás sumando presencia (1) y ausencia (0). En ese caso lo que obtienes es algo que se llama Modelo de Probabilidad Lineal.
Desde el punto de vista estadístico, donde la preocupación principal es la estimación de los errores estándar para hacer inferencia estadística este tipo de modelos no tiene mucho sentido y por eso se recomienda usar modelos logit.
Desde un punto de vista de ciencia de datos/machine learning/modelado predictivo la estimación de los errores estándar asociados a cada parámetro es mucho menos importante que el error promedio de clasificación en datos futuros (predicción). Entonces si un modelo de probabilidad lineal da buenas predicciones, adelante. Nada te impide que lo uses si en tu métrica de evaluación del modelo anda bien.
En el caso que presentas, sin embargo, esto no funcionaría, por lo dicho sobre las magnitudes implícitas en los números. Una alternativa sería un modelo logit multinomial. En R los implementa la librería nnet
y mlogit
.