La régression linéaire ne correspond pas bien

9

Je fais une régression linéaire en utilisant la fonction R lm:

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

entrez la description de l'image ici

mais ça ne va pas bien. Malheureusement, je ne peux pas comprendre le manuel.

Quelqu'un peut-il m'orienter dans la bonne direction pour mieux l'adapter?

En ajustant, je veux dire que je veux minimiser l'erreur quadratique moyenne (RMSE).


Edit : J'ai posté une question connexe (c'est le même problème) ici: Puis-je diminuer davantage le RMSE en fonction de cette fonctionnalité?

et les données brutes ici:

http://tny.cz/c320180d

sauf que sur ce lien x est ce qu'on appelle des erreurs sur la page actuelle ici, et il y a moins d'échantillons (1000 vs 3000 dans le tracé de la page actuelle). Je voulais simplifier les choses dans l'autre question.

Timothée HENRY
la source
4
R lm fonctionne comme prévu, le problème vient de vos données, c'est-à-dire que la relation linéaire n'est pas appropriée dans ce cas.
mpiktas
2
Pourriez-vous tracer quelle ligne vous pensez que vous devriez obtenir et pourquoi vous pensez que votre ligne a un MSE plus petit? Je note le mensonge de votre y entre 0 et 1, il semble donc que la régression linéaire serait tout à fait inappropriée pour ces données. Quelles sont les valeurs?
Glen_b -Reinstate Monica
2
Si les valeurs y sont des probabilités, vous ne voulez pas du tout de régression OLS.
Peter Flom
3
(désolé pourrais poster ceci avant) Ce qui vous semble être "un meilleur ajustement" ci-dessous est de minimiser (approximativement) les sommes des carrés des distances orthogonales, pas les distances verticales 'votre intuition se trompe. Vous pouvez vérifier le MSE approximatif assez facilement! Si les valeurs y sont des probabilités, vous seriez mieux servi par un modèle qui ne
sort
2
Il se pourrait que cette régression souffre de la présence de quelques valeurs aberrantes. Cela pourrait être un cas pour une régression robuste. en.wikipedia.org/wiki/Robust_regression
Yves Daoust

Réponses:

18

L'une des solutions les plus simples reconnaît que les changements parmi les probabilités qui sont faibles (comme 0,1) ou dont les compléments sont petits (comme 0,9) sont généralement plus significatifs et méritent plus de poids que les changements parmi les probabilités intermédiaires (comme 0,5).

Par exemple, un changement de 0,1 à 0,2 (a) double la probabilité tandis que (b) ne change la probabilité complémentaire que de 1/9 (en la faisant passer de 1-0,1 = 0,9 à 1-0,2 à 0,8), tandis qu'un changement de 0,5 à 0,6 (a) n'augmente la probabilité que de 20% tandis que (b) ne diminue la probabilité complémentaire que de 20%. Dans de nombreuses applications, ce premier changement est, ou du moins devrait être, considéré comme presque deux fois plus important que le second.

Dans toute situation où il serait tout aussi significatif d'utiliser une probabilité (que quelque chose se produise) ou son complément (c'est-à-dire la probabilité que quelque chose ne se produise pas), nous devons respecter cette symétrie.

Ces deux idées - de respecter la symétrie entre les probabilités p et leurs compléments 1pet d'exprimer les changements relativement plutôt qu'absolument - suggèrent que lorsque l'on compare deux probabilitésp et p nous devrions suivre à la fois leurs ratios p/p et les ratios de leurs compléments (1p)/(1p). Lors du suivi des ratios, il est plus simple d'utiliser des logarithmes, qui convertissent les ratios en différences. Ergo, un bon moyen d'exprimer une probabilitép à cet effet est d'utiliser

z=logplog(1p),
qui est connu comme les cotes de log ou logit dep. Cotes des bûches ajustéesz peut toujours être reconverti en probabilités en inversant le logit;
p=exp(z)/(1+exp(z)).
La dernière ligne du code ci-dessous montre comment cela se fait.

Ce raisonnement est assez général: il conduit à une bonne procédure initiale par défaut pour explorer tout ensemble de données impliquant des probabilités. (Il existe de meilleures méthodes disponibles, comme la régression de Poisson, lorsque les probabilités sont basées sur l'observation de ratios de «succès» par rapport au nombre d '«essais», car les probabilités basées sur plus d'essais ont été mesurées de manière plus fiable. Cela ne semble pas être le ici, où les probabilités sont basées sur des informations obtenues. On pourrait approcher l'approche de régression de Poisson en utilisant les moindres carrés pondérés dans l'exemple ci-dessous pour permettre des données plus ou moins fiables.)

Regardons un exemple.

Les figures

Le nuage de points à gauche montre un ensemble de données (similaire à celui de la question) tracé en termes de cotes logarithmiques. La ligne rouge correspond à ses moindres carrés ordinaires. Il a un faibleR2, indiquant beaucoup de dispersion et une forte "régression vers la moyenne": la droite de régression a une pente plus petite que l'axe principal de ce nuage de points elliptique. Il s'agit d'un cadre familier; il est facile à interpréter et à analyser en utilisant Rla lmfonction de ou l'équivalent.

Le nuage de points à droite exprime les données en termes de probabilités, telles qu'elles ont été enregistrées à l'origine. Le même ajustement est tracé: il semble maintenant incurvé en raison de la manière non linéaire dont les cotes de log sont converties en probabilités.

Dans le sens de l'erreur quadratique moyenne en termes de cotes logarithmiques, cette courbe est la mieux adaptée.

Par ailleurs, la forme approximativement elliptique du nuage à gauche et la façon dont il suit la ligne des moindres carrés suggèrent que le modèle de régression des moindres carrés est raisonnable: les données peuvent être décrites de manière adéquate par une relation linéaire - à condition que les cotes logarithmiques soient utilisées - et la variation verticale autour de la ligne est à peu près de la même taille indépendamment de l'emplacement horizontal (homoscédasticité). (Il y a quelques valeurs anormalement basses au milieu qui pourraient mériter un examen plus approfondi.) Évaluez cela plus en détail en suivant le code ci-dessous avec la commande plot(fit)pour voir des diagnostics standard. Cela est à lui seul une bonne raison d'utiliser les cotes de log pour analyser ces données au lieu des probabilités.


#
# 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)
whuber
la source
Merci beaucoup pour la réponse. J'aurai besoin de temps pour l'essayer.
Timothée HENRY du
Je rencontre une erreur en essayant votre code avec mes données, en essayant d'ajuster les cotes du journal: "Erreur dans lm.fit (x, y, offset = offset, singular.ok = singular.ok, ...): NA / NaN / Inf dans l'appel de fonction étrangère (arg 4) ".
Timothée HENRY du
Veuillez lire les commentaires dans le code: ils expliquent quel est le problème et comment y remédier.
whuber
6

Étant donné le biais dans les données avec x, la première chose évidente à faire est d'utiliser une régression logistique ( lien wiki ). Je suis donc avec whuber à ce sujet. Je dirai quexà lui seul montrera une forte signification mais n'expliquera pas la majeure partie de la déviance (l'équivalent de la somme totale des carrés dans un OLS). On pourrait donc suggérer qu'il existe une autre covariable en dehors dex qui facilite le pouvoir explicatif (par exemple, les personnes qui font la classification ou la méthode utilisée), Votre yles données sont déjà [0,1] cependant: savez-vous si elles représentent des probabilités ou des taux d'occurrence? Si c'est le cas, vous devriez essayer une régression logistique en utilisant votre non transforméy (avant qu'ils ne soient des ratios / probabilités).

L'observation de Peter Flom n'a de sens que si votre y n'est pas une probabilité. Vérifiez plot(density(y));rug(y)à différents seaux dexet voyez si vous voyez une distribution Bêta changeante ou exécutez simplement betareg. Notez que la distribution bêta est également une distribution de famille exponentielle et il devrait donc être possible de la modéliser avec glmdans R.

Pour vous donner une idée de ce que je voulais dire par régression logistique:

# 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)

Une régression logistique où le vrai modèle est $ log (\ frac {p} {1-p}) = 2-0,5x

EDIT: après avoir lu les commentaires:

Étant donné que «les valeurs y sont des probabilités d'appartenir à une certaine classe, obtenues à partir de classifications moyennes effectuées manuellement par des personnes», je recommande fortement de faire une régression logistique sur vos données de base. Voici un exemple:

Supposons que vous envisagez la probabilité qu'une personne accepte une proposition (y=1 se mettre d'accord, y=0 en désaccord) étant donné une incitation xentre 0 et 10 (pourrait être transformé en log, par exemple la rémunération). Deux personnes proposent l'offre aux candidats ("Jill and Jack"). Le vrai modèle est que les candidats ont un taux d'acceptation de base et qu'il augmente à mesure que la motivation augmente. Mais cela dépend aussi de qui propose l'offre (dans ce cas, nous disons que Jill a une meilleure chance que Jack). Supposons que, combinés, ils demandent 1000 candidats et collectent leurs données d'acceptation (1) ou de rejet (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) 

Du résumé, vous pouvez voir que le modèle s'adapte assez bien. La déviance estχn32 (std de χ2 est 2.df). Ce qui correspond et il bat un modèle avec une probabilité fixe (la différence d'écarts est de plusieurs centaines avecχ22). C'est un peu plus difficile à dessiner étant donné qu'il y a deux covariables ici mais vous avez l'idée.

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")

Jill dans Red Jack en bleu

Comme vous pouvez le voir, Jill a plus de facilité à obtenir un bon taux de réussite que Jack, mais cela disparaît à mesure que l'incitation augmente.

Vous devez essentiellement appliquer ce type de modèle à vos données d'origine. Si votre sortie est binaire, conservez le 1/0 s'il est multinomial vous avez besoin d'une régression logistique multinomiale. Si vous pensez que la source supplémentaire de variance n'est pas le collecteur de données, ajoutez un autre facteur (ou variable continue) selon ce que vous jugez logique pour vos données. Les données viennent en premier, deuxième et troisième, seulement alors le modèle entre en jeu.

Hans Roggeman
la source
Un commentaire du PO, "Les valeurs y sont des probabilités d'appartenir à une certaine classe, obtenues à partir de classifications moyennes effectuées manuellement par des personnes", suggère que la régression logistique serait inappropriée pour ces données - même si cela pourrait être une excellente solution pour le des données brutes (comme suggéré dans votre premier paragraphe), selon ce que sont les "classifications" et comment la "moyenne" s'est produite. Appliqué aux données présentées dans la question, glmproduit une ligne incurvée relativement plate qui ressemble remarquablement à la ligne indiquée dans la question.
whuber
Je vous remercie. Et oui, y est une probabilité. J'ai également publié les données brutes dans une question connexe: stats.stackexchange.com/questions/83576/… , bien que j'aie appelé x ce que j'ai appelé log (x) dans l'autre question ...
Timothée HENRY
J'aurais aimé le savoir avant d'acquérir un échantillon de votre image, LOL!
whuber
5

Le modèle de régression linéaire n'est pas bien adapté aux données. On pourrait s'attendre à tirer quelque chose comme ce qui suit de la régression:

entrez la description de l'image ici

mais en réalisant ce que fait OLS, il est évident que ce n'est pas ce que vous obtiendrez. Une interprétation graphique des moindres carrés ordinaires est qu'elle minimise la distance verticale au carré entre la ligne (hyperplan) et vos données. De toute évidence, la ligne violette que j'ai superposée a d'énormes résidus dex(7,4.5) et encore de l'autre côté de 3. C'est pourquoi la ligne bleue est mieux adaptée que la violette.

pkofod
la source
@pkofod Oui, je vois. J'ai donc supprimé mon commentaire (je savais que vous saviez qu'il était carré; mais d'autres lecteurs pourraient ne pas le faire).
Peter Flom
1
La régression censurée est différente de la régression avec une variable dépendante limitée à une plage connue fixe. Ces données ne sont pas censurées et la régression censurée ne fera rien de différent avec elles que la régression ordinaire.
whuber
Ouais, ma mauvaise. Supprimé cette partie.
pkofod
4

Puisque Y est borné par 0 et 1, la régression des moindres carrés ordinaires n'est pas bien adaptée. Vous pouvez essayer la régression bêta. Il Ry a le betaregpaquet.

Essayez quelque chose comme ça

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

Plus d'informations

EDIT: Si vous souhaitez un compte rendu complet de la régression bêta, ses avantages et ses inconvénients, voir Un meilleur presse-citron: régression de vraisemblance maximale avec des variables dépendantes distribuées bêta par Smithson et Verkuilen

Peter Flom
la source
4
Quel modèle est en betaregtrain de mettre en œuvre? Quelles sont ses hypothèses et pourquoi est-il raisonnable de supposer qu'elles s'appliquent à ces données?
whuber
2
@whuber C'est une bonne question! Le modèle est défini aux pages 3 et 4 de cette vignette . Il est basé sur une densité bêta reparamétrisée en termes de paramètres de moyenne et de précision (qui peuvent tous deux être modélisés, chacun avec sa propre fonction de lien), et un ensemble de fonctions de lien identiques à celles utilisées pour les modèles binomiaux (et un de plus). Il est monté par ML et fonctionne très comme un GLM.
Glen_b -Reinstate Monica
2
@whuber Les modèles bêta conditionnels sont courants pour les données de composition et d'autres proportions ou probabilités de type continu. Je ne sais pas si les hypothèses pour de tels modèles conviennent à ces données (je ne sais pas quelles sont les données, ce qui serait ma première préoccupation avant de suggérer un modèle moi-même), mais même juste à partir de l'intrigue, j'imagine que ils conviendraient ainsi que d'autres modèles suggérés ici. Il y a un certain nombre de modèles dans les réponses qui ne semblent pas être mieux justifiés que la suggestion de Peter, certains avec des hypothèses (pas toujours énoncées) qui semblent plus difficiles à justifier.
Glen_b -Reinstate Monica
1
Merci, @Glen_b. Je ne conteste pas la suggestion de Peter - j'essaie seulement de la comprendre, car je n'ai jamais utilisé de régression bêta auparavant, et j'imagine que de nombreux futurs lecteurs seraient dans la même situation. (Je connais suffisamment les autres modèles mentionnés dans ce fil pour comprendre leurs hypothèses et leurs éventuelles lacunes!) Il serait donc agréable de voir cette réponse inclure au moins un bref compte rendu des hypothèses et des raisons de recommander cette solution.
whuber
1
Ah, oui, je me suis lié à ce document dans une réponse moi-même à un moment donné. Smithson (l'un des auteurs) a le document sur sa page Web . Le matériel supplémentaire est lié ici .
Glen_b -Reinstate Monica
1

Vous voudrez peut-être d'abord savoir exactement ce que fait un modèle linéaire. Il essaie de modéliser une relation de la forme

Yi=a+bXi+ϵi

Où le ϵisatisfaire à certaines conditions (hétéroskédasticité, variance uniforme et indépendance - wikipedia est un bon début si cela ne vous dit rien). Mais même si ces conditions sont vérifiées, il n'y a absolument aucune garantie que ce sera le "meilleur ajustement" dans le sens que vous recherchez: OLS essaie simplement de minimiser les erreurs dans la direction Y, ce qu'il fait dans votre mais ce n'est pas ce qui semble être le mieux adapté.

Si un modèle linéaire est vraiment ce que vous recherchez, vous pouvez essayer de transformer un peu vos variables afin que OLS puisse en effet être ajusté, ou simplement essayer un autre modèle. Vous voudrez peut-être examiner PCA ou CCA, ou si vous êtes vraiment déterminé à utiliser un modèle linéaire, essayez la solution des moindres carrés , qui pourrait donner un meilleur "ajustement", car elle permet des erreurs dans les deux sens.

Youloush
la source
Je pensais que lm cherchait un minimum de "moindres carrés totaux" pour une fonction linéaire (a + b * x + epsilon). Je suis perdu.
Timothée HENRY du
1
lm, comme vous l'avez utilisé, minimise la somme des résidus au carré, c'est-à-dire (yabx)2pour chaque point de données, appelé OLS (moindres carrés ordinaires). Je n'ai pas pu trouver une belle image pour OLS linéaire, mais peut - être que celle-ci est encore illustrative pour vous. OLS minimise le carré des lignes vertes, lm le fait avec une ligne. En comparaison, regardez un linéaire des moindres carrés totaux .
Roland