[R] Predicting probabilities from a logistic regression by hand (in code)
erikpukinskis
erik at snowedin.net
Fri Jan 17 03:53:56 CET 2014
Thanks for looking at this, I've been tearing my hair out for a day or so
now.
I have done a multiple variable logistic regression in R, and obtained my
coefficients. I am able to make predictions for the training data in R
without problem. But now I would like to create a prediction model in Ruby
(that was the original point of doing the regression) and I'm having some
trouble.
Basically, my equation is:
predicted_logit = K + v1*c1 + v2*c2 + ... vn*cn
odds_ratio = e^predicted_logit/(1+e^predicted_logit)
But it always seems to either give 1.0 or 0.0! The output of predict() in R
is generally something nice and soft like 0.5578460!
I realize not everyone knows Ruby, but I'll include my code here for
reference:
# These are the coefficients that R gives me from my logistic regression:
intercept = 0.2700309
coefficients = {
high: 1.0136028,
low: 1.0016712,
germ_mean: 1.0233327,
gdds: 0.9990283,
early_gdds: 0.9986464,
mid_gdds: 1.0002979,
late_gdds: 0
}
# And this is what R predicts for one datum:
#
# outcome high low germ_mean gdds early_gdds mid_gdds late_gdds p_success
# 1 1 73 28 40 119 0 91 28 0.5578460
# ...
# So to get my own p_success, first I multiply each coefficient by it's
input data
period = {:high=>73, :low=>28, :germ_mean=>40, :gdds=>119, :early_gdds=>0,
:mid_gdds=>91, :late_gdds=>28}
products = coefficients.map {|name,value| period[name]*value }
# Then I add those together and add that to the intercept
predicted_logit = intercept + products.sum
# Then my probability should be e^predicted_logit over 1 +
e^predicted_logit:
odds_ratio = Math.exp(predicted_logit) / (1 + Math.exp(predicted_logit))
# But the odds ratio comes out as 1.0, not 0.5578460 like R predicts.
--
View this message in context: http://r.789695.n4.nabble.com/Predicting-probabilities-from-a-logistic-regression-by-hand-in-code-tp4683713.html
Sent from the R help mailing list archive at Nabble.com.
More information about the R-help
mailing list