La regresión lineal no se ajusta bien

Question

Hago una regresión lineal utilizando la función lm de R:

x = log(errors)
plot(x,y)
lm.result = lm(formula = y ~ x)
abline(lm.result, col="blue") # showing the "fit" in blue

pero no encaja bien. Lamentablemente no puedo entender el manual.

¿Puede alguien indicarme la dirección correcta para encajar mejor esto?

Por ajuste quiero decir que quiero minimizar el error medio cuadrático (RMSE).

---

Editar : He publicado una pregunta relacionada (es el mismo problema) aquí: ¿Puedo disminuir aún más el RMSE basándome en esta característica?

y los datos en bruto aquí:

http://tny.cz/c320180d

excepto que en ese enlace x es lo que se llama errores en la presente página aquí, y hay menos muestras (1000 frente a 3000 en el gráfico de la presente página). Quería simplificar las cosas en la otra pregunta.

Answer 1

Una de las soluciones más sencillas reconoce que los cambios entre las probabilidades que son pequeñas (como 0,1) o cuyos complementos son pequeños (como 0,9) suelen ser más significativos y merecen más peso que los cambios entre las probabilidades medias (como 0,5).

Por ejemplo, un cambio de 0,1 a 0,2 (a) duplica la probabilidad mientras que (b) cambia la probabilidad complementaria sólo en 1/9 (bajando de 1-0,1 = 0,9 a 1-0,2 a 0,8), mientras que un cambio de 0,5 a 0,6 (a) aumenta la probabilidad sólo en un 20% mientras que (b) disminuye la probabilidad complementaria sólo en un 20%. En muchas aplicaciones, el primer cambio se considera, o al menos debería considerarse, casi el doble que el segundo.

En cualquier situación en la que tenga el mismo sentido utilizar una probabilidad (de que algo ocurra) o su complemento (es decir, la probabilidad de que ese algo no ocurra), debemos respetar esta simetría.

Estas dos ideas -de respetar la simetría entre las probabilidades $p$ y sus complementos $1-p$ y de expresar los cambios de forma relativa y no absoluta, sugieren que al comparar dos probabilidades $p$ y $p'$ deberíamos seguir sus dos ratios $p'/p$ y las relaciones de sus complementos $(1-p)/(1-p')$ . Para el seguimiento de las relaciones es más sencillo utilizar los logaritmos, que convierten las relaciones en diferencias. Ergo, una buena forma de expresar una probabilidad $p$ para este fin es utilizar $$z = \log p - \log(1-p),$$ que se conoce como el probabilidades de registro o logit de $p$ . Probabilidades logarítmicas ajustadas $z$ siempre se pueden volver a convertir en probabilidades invirtiendo el logit; $$p = \exp(z)/(1+\exp(z)).$$ La última línea del código siguiente muestra cómo se hace esto.

Este razonamiento es bastante general: conduce a un buen procedimiento inicial por defecto para explorar cualquier conjunto de datos que impliquen probabilidades. (Existen mejores métodos, como la regresión de Poisson, cuando las probabilidades se basan en la observación de las proporciones de "éxitos" con respecto al número de "ensayos", porque las probabilidades basadas en más ensayos se han medido con mayor fiabilidad. Este no parece ser el caso aquí, donde las probabilidades se basan en la información obtenida. Se podría aproximar el enfoque de la regresión de Poisson utilizando los mínimos cuadrados ponderados en el ejemplo siguiente para tener en cuenta los datos que son más o menos fiables).

Veamos un ejemplo.

El gráfico de dispersión de la izquierda muestra un conjunto de datos (similar al de la pregunta) representado en términos de probabilidades logarítmicas. La línea roja es su ajuste por mínimos cuadrados ordinarios. Tiene una baja $R^2$ que indica mucha dispersión y una fuerte "regresión a la media": la línea de regresión tiene una pendiente menor que el eje mayor de esta nube de puntos elíptica. Se trata de un escenario familiar; es fácil de interpretar y analizar utilizando R 's lm o su equivalente.

El gráfico de dispersión de la derecha expresa los datos en términos de probabilidades, tal y como se registraron originalmente. El mismo Ahora se ve curvado debido a la forma no lineal en que las probabilidades logarítmicas se convierten en probabilidades.

En el sentido del error medio cuadrático en términos de probabilidades logarítmicas, esta curva es la mejor encaja.

Por cierto, la forma aproximadamente elíptica de la nube de la izquierda y la forma en que sigue la línea de mínimos cuadrados sugieren que el modelo de regresión de mínimos cuadrados es razonable: los datos pueden describirse adecuadamente mediante una relación lineal proporcionó y la variación vertical en torno a la línea es aproximadamente del mismo tamaño, independientemente de la ubicación horizontal (homocedasticidad). (Hay algunos valores inusualmente bajos en el centro que podrían merecer un análisis más detallado). Evalúe esto con más detalle siguiendo el código siguiente con el comando plot(fit) para ver algunos diagnósticos estándar. Esto por sí solo es una razón de peso para utilizar las probabilidades logarítmicas para analizar estos datos en lugar de las probabilidades.

---

#
# Read the data from a table of (X,Y) = (X, probability) pairs.
#
x <- read.table("F:/temp/data.csv", sep=",", col.names=c("X", "Y"))
#
# Define functions to convert between probabilities `p` and log odds `z`.
# (When some probabilities actually equal 0 or 1, a tiny adjustment--given by a positive
# value of `e`--needs to be applied to avoid infinite log odds.)
#
logit <- function(p, e=0) {x <- (p-1/2)*(1-e) + 1/2; log(x) - log(1-x)}
logistic <- function(z, e=0) {y <- exp(z)/(1 + exp(z)); (y-1/2)/(1-e) + 1/2}
#
# Fit the log odds using least squares.
#
b <- coef(fit <- lm(logit(x$Y) ~ x$X))
#
# Plot the results in two ways.
#
par(mfrow=c(1,2))
plot(x$X, logit(x$Y), cex=0.5, col="Gray",
     main="Least Squares Fit", xlab="X", ylab="Log odds")
abline(b, col="Red", lwd=2)

plot(x$X, x$Y, cex=0.5, col="Gray",
     main="LS Fit Re-expressed", xlab="X", ylab="Probability")
curve(logistic(b[1] + b[2]*x), col="Red", lwd=2, add=TRUE)

Answer 2

Dado el sesgo de los datos con x, lo primero que hay que hacer es utilizar una regresión logística ( enlace wiki ). Así que estoy con Whuber en esto. Diré que $x$ por sí sola mostrará una fuerte significación pero no explicará la mayor parte de la desviación (el equivalente a la suma total de cuadrados en un MCO). Así que se podría sugerir que hay otra covariable aparte de $x$ que ayuda al poder explicativo (por ejemplo, las personas que realizan la clasificación o el método utilizado), Su $y$ los datos ya son [0,1]: ¿sabe si representan probabilidades o ratios de ocurrencia? Si es así, debería intentar una regresión logística utilizando sus datos no transformados $y$ (antes de que sean ratios/probabilidades).

La observación de Peter Flom sólo tiene sentido si su y no es una probabilidad. Compruebe plot(density(y));rug(y) en diferentes cubos de $x$ y vea si ve una distribución Beta cambiante o simplemente ejecute betareg . Obsérvese que la distribución beta es también una distribución de la familia exponencial y, por tanto, debería ser posible modelarla con glm en R.

Para que te hagas una idea de lo que quiero decir con la regresión logística:

# the 'real' relationship where y is interpreted as the probability of success
y = runif(400)
x = -2*(log(y/(1-y)) - 2) + rnorm(400,sd=2) 
glm.logit=glm(y~x,family=binomial); summary(glm.logit) 
plot(y ~ x); require(faraway); grid()
points(x,ilogit(coef(glm.logit) %*% rbind(1.0,x)),col="red")
tt=runif(400)  # an example of your untransformed regression
newy = ifelse(tt < y, 1, 0)
glm.logit=glm(newy~x,family=binomial); summary(glm.logit) 

# if there is not a good match in your tail probabilities try different link function or oversampling with correction (will be worse here, but perhaps not in your data)
glm.probit=glm(y~x,family=binomial(link=probit)); summary(glm.probit)
glm.cloglog=glm(y~x,family=binomial(link=cloglog)); summary(glm.cloglog)

EDIT: después de leer los comentarios:

Dado que "Los valores de y son probabilidades de pertenecer a una determinada clase, obtenidas de promediar clasificaciones hechas manualmente por personas", recomiendo encarecidamente hacer una regresión logística sobre tus datos base. He aquí un ejemplo:

Supongamos que se trata de la probabilidad de que alguien acepte una propuesta ( $y=1$ estar de acuerdo, $y=0$ no están de acuerdo) dado un incentivo $x$ entre 0 y 10 (podría transformarse en logaritmo, por ejemplo, la remuneración). Hay dos personas que proponen la oferta a los candidatos ("Jill y Jack"). El modelo real es que los candidatos tienen una tasa de aceptación base y que ésta aumenta a medida que aumenta el incentivo. Pero también depende de quién proponga la oferta (en este caso decimos que Jill tiene más posibilidades que Jack). Supongamos que en conjunto preguntan a 1000 candidatos y recogen sus datos de aceptación (1) o rechazo (0).

require(faraway)
people = c("Jill","Jack")
proposer = sample(people,1000,replace=T)
incentive = runif(1000, min = 0, max =10)
noise = rnorm(1000,sd=2)
# base probability of agreeing is about 12% (ilogit(-2))
agrees = ilogit(-2 + 1*incentive + ifelse(proposer == "Jill", 0 , -0.75) + noise) 
tt = runif(1000)
observedAgrees = ifelse(tt < agrees,1,0)
glm.logit=glm(observedAgrees~incentive+proposer,family=binomial); summary(glm.logit) 

En el resumen se puede ver que el modelo se ajusta bastante bien. La desviación es $\chi^2_{n-3}$ (std de $\chi^2$ es $\sqrt{2.df}$ ). Que se ajusta y supera a un modelo con una probabilidad fija (la diferencia de desviaciones es de varios cientos con $\chi^2_{2}$ ). Es un poco más difícil de dibujar dado que hay dos covariables aquí, pero se entiende la idea.

xs = coef(glm.logit) %*% rbind(1,incentive,as.factor(proposer))
ys = as.vector(unlist(ilogit(xs)))
plot(ys~ incentive, type="n"); require(faraway); grid()
points(incentive[proposer == "Jill"],ys[proposer == "Jill"],col="red")
points(incentive[proposer == "Jack"],ys[proposer == "Jack"],col="blue")

Como puedes ver, a Jill le resulta más fácil conseguir un buen porcentaje de aciertos que a Jack, pero eso desaparece a medida que aumenta el incentivo.

Básicamente, debe aplicar este tipo de modelo a sus datos originales. Si tu salida es binaria, mantén el 1/0 si es multinomial necesitas una regresión logística multinomial. Si crees que la fuente extra de varianza no es el colector de datos, añade otro factor (o variable continua) lo que creas que tiene sentido para tus datos. Los datos son lo primero, lo segundo y lo tercero, sólo entonces entra en juego el modelo.

Answer 3

El modelo de regresión lineal no se adapta bien a los datos. Se podría esperar obtener algo como lo siguiente de la regresión:

pero al darse cuenta de lo que hace OLS, es obvio que esto no es lo que se obtiene. Una interpretación gráfica de los mínimos cuadrados ordinarios es que minimiza el cuadrado vertical distancia entre la línea (hiperplano) y sus datos. Obviamente, la línea púrpura que he superpuesto tiene unos residuos enormes de $x\in (-7,4.5)$ y de nuevo al otro lado del 3. Por eso la línea azul se ajusta mejor que la púrpura.

Answer 4

Dado que Y está delimitado por 0 y 1, la regresión por mínimos cuadrados ordinarios no es adecuada. Puedes probar con la regresión beta. En R está el betareg paquete.

Intenta algo como esto

install.packages("betareg")
library(betareg)
betamod1 <- betareg(y~x, data = DATASETNAME)

más información

EDIT: Si desea una explicación completa de la regresión beta, sus ventajas y desventajas, consulte Un mejor exprimidor de limones: Regresión de máxima verosimilitud con variables dependientes con distribución beta por Smithson y Verkuilen

Answer 5

Es posible que primero quiera saber qué hace exactamente un modelo lineal. Intenta modelar una relación de la forma

$$Y_i = a + bX_i + \epsilon_i$$

Donde el $\epsilon_i$ satisfacen ciertas condiciones (heteroscedasticidad, varianza uniforme e independencia - la wikipedia es un buen comienzo si no le suena). Pero incluso si se verifican estas condiciones, no hay absolutamente ninguna garantía de que esto será un "mejor ajuste" en el sentido que usted está buscando : OLS sólo trata de minimizar el error en la dirección Y, que es lo que está haciendo en su caso, pero que no es lo que parece el mejor ajuste.

Si lo que realmente busca es un modelo lineal, puede intentar transformar un poco sus variables para que se puedan ajustar las MCO, o simplemente probar otro modelo. Puede que quiera buscar en PCA o CCA, o si está realmente empeñado en usar un modelo lineal, pruebe el mínimos cuadrados totales solución, que podría dar un mejor "ajuste", ya que permite errores en ambas direcciones.