如何获取维基百科R中逻辑回归模仿的ggplot图?
How to obtain the ggplot graph for the logistic regression imitation in Wikipedia's example in R?
(添加了可重现的示例)
我试图模仿维基百科的 "Probability of passing an exam versus hours of study" 逻辑回归示例 here:
我无法在那个页面中获得相同的 ggplot 图,我也不知道为什么。
df <- data.frame(hour=c(0.50,0.75,1.00,1.25,1.50,1.75,1.75,2.00,2.25,2.50,2.75,3.00,3.25,3.50,4.00,4.25,4.50,4.75,5.00,5.50), pass=c(0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,1,1,1,1,1))
df
# hour pass
#1 0.50 0
#2 0.75 0
#3 1.00 0
#4 1.25 0
#5 1.50 0
#6 1.75 0
#7 1.75 1
#8 2.00 0
#9 2.25 1
#10 2.50 0
#11 2.75 1
#12 3.00 0
#13 3.25 1
#14 3.50 0
#15 4.00 1
#16 4.25 1
#17 4.50 1
#18 4.75 1
#19 5.00 1
#20 5.50 1
df$pass <- as.factor(df$pass)
my_fit <- glm(df$pass ~ df$hour, data=df, na.action=na.exclude, family="binomial")
summary(my_fit)
非 GGPLOT 情节完美运行:
my_table <- summary(my_fit)
my_table$coefficients[,1] <- invlogit(coef(my_fit))
my_table
anova(my_fit)
library(pscl); pR2(my_fit) # for McFadden rho^2
plot(df$hour, df$pass, xlab="x", ylab="logit values")
LinearPredictions <- predict(my_fit); LinearPredictions
# LinearPredictions is NOT equal to 0.01666 + 0.81827*(1:20)
# LinearPredictions is NOT equal to -4.0777+1.5046*(1:20)
# LinearPredictions are equal to what (I couldn't solve)?
EstimatedProbability.hat <- exp(LinearPredictions)/(1 + exp(LinearPredictions))
EstimatedProbability.hat
EstimatedProbability <- c(0.25, 0.50, 0.75) # Estimated probabilities for which their x levels are wanted to be found
HoursStudied <- (log(EstimatedProbability/(1- EstimatedProbability)) - my_fit$coefficients[1])/ my_fit$coefficients[2]
HoursStudied.summary <- data.frame(EstimatedProbability, HoursStudied)
HoursStudied.summary
plot(df$hour, EstimatedProbability.hat, xlab="studying hours", ylab="estimated probability (pass)") # , xlim=c(0,6), ylim=c(0,1)
# Add red curve
lines(df$hour, EstimatedProbability.hat, lty=1, col="red")
# Vertical dashes
segments(x0=HoursStudied.summary$HoursStudied, y0=0, x1=HoursStudied.summary$HoursStudied, y1=HoursStudied.summary$EstimatedProbability,
lty=2, col=c("darkblue","darkred","darkgreen"))
# Horizontal dashes
segments(x0=0, y0=HoursStudied.summary$EstimatedProbability, x1=HoursStudied.summary$HoursStudied,
y1=HoursStudied.summary$EstimatedProbability, lty=2, col=c("darkblue","darkred","darkgreen"))
legend("bottomright", legend=c("HS0.25", "HS0.50", "HS0.75"), lty=2, col=c("darkblue","darkred","darkgreen"), bty="n", cex=0.75)
GGPLOT 绘图失败:
我试图在 ggplot
中做同样的事情,但失败了:
df$EstimatedProbabilities <- EstimatedProbability.hat; df
HoursStudied.summary$group <- c('HS0.25','HS0.50','HS0.75')
library(ggplot2)
ggplot(df, aes(x=hour, y=df$pass)) +
geom_point() +
geom_line(aes(y=EstimatedProbabilities), colour="black") +
geom_segment(data=HoursStudied.summary, aes(y=EstimatedProbability,
xend=HoursStudied, yend=EstimatedProbability, col=group), x=-Inf, linetype="dashed") +
geom_segment(data=HoursStudied.summary, aes(x=HoursStudied,
xend=HoursStudied, yend=EstimatedProbability, col=group), y=-Inf, linetype="dashed")
问题:ggplot
曲线与plot
曲线相同,但是它把整个函数绘制在y=0线下方。为什么?
编辑: 您需要 df$pass
是数字,而不是因数。我也不会在最初的 ggplot
调用中映射任何美学,而只是在 geom_point
和 geom_line
调用中传递它们。
df$pass <- as.numeric(df$pass) - 1
ggplot(df) +
geom_point(aes(x=hour,y=pass)) +
geom_line(aes(x=hour,y=EstimatedProbabilities)) +
geom_segment(data=HoursStudied.summary, aes(y=EstimatedProbability, xend=HoursStudied, yend=EstimatedProbability, col=group), x=-Inf, linetype="dashed") +
geom_segment(data=HoursStudied.summary, aes(x=HoursStudied, xend=HoursStudied, yend=EstimatedProbability, col=group), y=-Inf, linetype="dashed")
这个问题使 geom_smooth
可以变得简单的事情变得复杂。请注意,预测是 type = "response"
,在 this post to CrossValidated.
之后
my_fit <- glm(pass ~ hour, data = df, na.action = na.exclude,
family = "binomial")
pred <- predict(my_fit, type = "response")
pred_df <- data.frame(hour = df$hour, pred)
library(ggplot2)
ggplot(df, aes(x = hour, y = pass)) +
geom_point() +
geom_smooth(method = "glm",
method.args = list(family = "binomial"),
se = FALSE) +
geom_point(data = pred_df, aes(x = hour, y = pred), colour = "blue") +
geom_hline(data = data.frame(c(0.25, 0.50, 0.75)),
aes(yintercept = c(0.25, 0.50, 0.75)),
colour = "darkgrey", linetype = "dashed")
(添加了可重现的示例)
我试图模仿维基百科的 "Probability of passing an exam versus hours of study" 逻辑回归示例 here:
我无法在那个页面中获得相同的 ggplot 图,我也不知道为什么。
df <- data.frame(hour=c(0.50,0.75,1.00,1.25,1.50,1.75,1.75,2.00,2.25,2.50,2.75,3.00,3.25,3.50,4.00,4.25,4.50,4.75,5.00,5.50), pass=c(0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,1,1,1,1,1))
df
# hour pass
#1 0.50 0
#2 0.75 0
#3 1.00 0
#4 1.25 0
#5 1.50 0
#6 1.75 0
#7 1.75 1
#8 2.00 0
#9 2.25 1
#10 2.50 0
#11 2.75 1
#12 3.00 0
#13 3.25 1
#14 3.50 0
#15 4.00 1
#16 4.25 1
#17 4.50 1
#18 4.75 1
#19 5.00 1
#20 5.50 1
df$pass <- as.factor(df$pass)
my_fit <- glm(df$pass ~ df$hour, data=df, na.action=na.exclude, family="binomial")
summary(my_fit)
非 GGPLOT 情节完美运行:
my_table <- summary(my_fit)
my_table$coefficients[,1] <- invlogit(coef(my_fit))
my_table
anova(my_fit)
library(pscl); pR2(my_fit) # for McFadden rho^2
plot(df$hour, df$pass, xlab="x", ylab="logit values")
LinearPredictions <- predict(my_fit); LinearPredictions
# LinearPredictions is NOT equal to 0.01666 + 0.81827*(1:20)
# LinearPredictions is NOT equal to -4.0777+1.5046*(1:20)
# LinearPredictions are equal to what (I couldn't solve)?
EstimatedProbability.hat <- exp(LinearPredictions)/(1 + exp(LinearPredictions))
EstimatedProbability.hat
EstimatedProbability <- c(0.25, 0.50, 0.75) # Estimated probabilities for which their x levels are wanted to be found
HoursStudied <- (log(EstimatedProbability/(1- EstimatedProbability)) - my_fit$coefficients[1])/ my_fit$coefficients[2]
HoursStudied.summary <- data.frame(EstimatedProbability, HoursStudied)
HoursStudied.summary
plot(df$hour, EstimatedProbability.hat, xlab="studying hours", ylab="estimated probability (pass)") # , xlim=c(0,6), ylim=c(0,1)
# Add red curve
lines(df$hour, EstimatedProbability.hat, lty=1, col="red")
# Vertical dashes
segments(x0=HoursStudied.summary$HoursStudied, y0=0, x1=HoursStudied.summary$HoursStudied, y1=HoursStudied.summary$EstimatedProbability,
lty=2, col=c("darkblue","darkred","darkgreen"))
# Horizontal dashes
segments(x0=0, y0=HoursStudied.summary$EstimatedProbability, x1=HoursStudied.summary$HoursStudied,
y1=HoursStudied.summary$EstimatedProbability, lty=2, col=c("darkblue","darkred","darkgreen"))
legend("bottomright", legend=c("HS0.25", "HS0.50", "HS0.75"), lty=2, col=c("darkblue","darkred","darkgreen"), bty="n", cex=0.75)
GGPLOT 绘图失败:
我试图在 ggplot
中做同样的事情,但失败了:
df$EstimatedProbabilities <- EstimatedProbability.hat; df
HoursStudied.summary$group <- c('HS0.25','HS0.50','HS0.75')
library(ggplot2)
ggplot(df, aes(x=hour, y=df$pass)) +
geom_point() +
geom_line(aes(y=EstimatedProbabilities), colour="black") +
geom_segment(data=HoursStudied.summary, aes(y=EstimatedProbability,
xend=HoursStudied, yend=EstimatedProbability, col=group), x=-Inf, linetype="dashed") +
geom_segment(data=HoursStudied.summary, aes(x=HoursStudied,
xend=HoursStudied, yend=EstimatedProbability, col=group), y=-Inf, linetype="dashed")
问题:ggplot
曲线与plot
曲线相同,但是它把整个函数绘制在y=0线下方。为什么?
编辑: 您需要 df$pass
是数字,而不是因数。我也不会在最初的 ggplot
调用中映射任何美学,而只是在 geom_point
和 geom_line
调用中传递它们。
df$pass <- as.numeric(df$pass) - 1
ggplot(df) +
geom_point(aes(x=hour,y=pass)) +
geom_line(aes(x=hour,y=EstimatedProbabilities)) +
geom_segment(data=HoursStudied.summary, aes(y=EstimatedProbability, xend=HoursStudied, yend=EstimatedProbability, col=group), x=-Inf, linetype="dashed") +
geom_segment(data=HoursStudied.summary, aes(x=HoursStudied, xend=HoursStudied, yend=EstimatedProbability, col=group), y=-Inf, linetype="dashed")
这个问题使 geom_smooth
可以变得简单的事情变得复杂。请注意,预测是 type = "response"
,在 this post to CrossValidated.
my_fit <- glm(pass ~ hour, data = df, na.action = na.exclude,
family = "binomial")
pred <- predict(my_fit, type = "response")
pred_df <- data.frame(hour = df$hour, pred)
library(ggplot2)
ggplot(df, aes(x = hour, y = pass)) +
geom_point() +
geom_smooth(method = "glm",
method.args = list(family = "binomial"),
se = FALSE) +
geom_point(data = pred_df, aes(x = hour, y = pred), colour = "blue") +
geom_hline(data = data.frame(c(0.25, 0.50, 0.75)),
aes(yintercept = c(0.25, 0.50, 0.75)),
colour = "darkgrey", linetype = "dashed")