为我的逻辑回归模型生成逻辑曲线
Producing logistic curve for my logistic regression model
我想编写代码来绘制我的逻辑回归模型,即 "S" 形状的逻辑曲线。我有 两个独立的协变量 ,请问如何做到这一点?我附上了我的数据集和我的模型的代码。先感谢您。
239 0.72 1
324.6 0.83 1
331.8 0.95 1
334.3 0.83 1
259.7 0.89 1
212.3 0.88 1
204.7 0.65 1
253.86 0.75 1
258.94 0.85 1
329.66 0.95 0
469.68 1.46 0
459.74 1.11 0
293.2 0.64 0
297.88 0.98 0
267.9 0.82 0
374.1 1.29 0
333.62 0.74 0
dat <- read.table("data.txt")
colnames(dat)<-c("press","v","gender")
# logostic regression
dat$gender <- factor(dat$gender)
mylogit<- glm(gender~press+v,data=dat,family="binomial")
summary(mylogit)
######## the code below are irrelevant to making plot, ignore if you want
mylogit$fitted.values
newdat <- data.frame(t(c(300,0.1)))
colnames(newdat)<-c("press","v")
# this is your new dataset, we name it as "newdat"
pred <- predict(mylogit,newdata = newdat,type="response")
pred # the probability of being in class 1 will stored in this object
pred <- predict(mylogit,newdata = dat,type="response")
pred # the probability of being in class 1 will stored in this object
# accuracy
dat$pred <- 0
factor(dat$pred)
dat$pred[which(pred>0.5)] <- 1
table(dat$gender,dat$pred)
您有 2 个连续的非分类变量,因此逻辑曲线将是 3D 曲线。我给大家提供两种展示方式。
- 使用
persp
函数生成真实的3D平滑曲线;
- 固定
v
为多个值,然后生成多个二维logistic曲线(你称之为"S"-shape曲线)。
3D曲线
press_grid <- seq(200, 480, by = 5)
v_grid <- seq(0.6, 1.5, by = 0.1)
newdat <- data.frame(press = rep(press_grid, times = length(v_grid)), v = rep(v_grid, each = length(press_grid)))
pred <- predict.glm(mylogit, newdata = newdat, type="response")
z <- matrix(pred, length(press_grid))
persp(press_grid, v_grid, z, xlab = "pressure", ylab = "velocity", zlab = "predicted probability", main = "logistic curve (3D)", theta = 30, phi = 20)
您需要先生成一个二维网格。 newdat
保存此网格,您可以执行 plot(newdat)
来查看此网格。然后通过调用 predict.glm(..., type = "response")
在此网格上进行预测。结果 pred
是一个向量。要绘制它,将其转换为矩阵 z
,然后调用 persp
进行 3D 绘图。 xlab
、ylab
、zlab
是三轴的标签。参数 theta
和 phi
用于调整视角。
在上面,press
和v
的边际网格是基于你的原始数据的范围:range(dat$press)
和range(dat$v)
。我们不做超出这个范围的预测。但即使在这个范围内,您也只有 17 个观测值。所以你仍然需要对情节持怀疑态度。
这是曲线:
二维曲线
这个玩具函数对于制作 2D 曲线很有用,v
固定为某个级别:
curve_2D_fix_v <- function(model, v = 1, press_grid = seq(200, 480, by = 5), add = FALSE, col = "black") {
newdat <- data.frame(press = press_grid, v = v)
pred <- predict.glm(model, newdat, type = "response")
if (add) lines(press_grid, pred, col = col) else {
plot(press_grid, pred, xlab = "pressure", ylab = "predicted probability", type = "l", col = col, main = "logistic curve (2D)")
abline(h = c(0, 0.5, 1), lty = 2, col = col)
}
}
如果add = FALSE
,则开启新的绘图window;虽然它是 TRUE
,但它绘制在之前的 window 上(但你有责任确保有这样一个 window!)二维图提供了更多信息,因为你可以添加一个0、0.5 和 1 处的水平线。
让我们开始吧:
curve_2D_fix_v(mylogit, v = 0.4, add = FALSE, col = "black")
curve_2D_fix_v(mylogit, v = 0.6, add = TRUE, col = "red")
curve_2D_fix_v(mylogit, v = 0.8, add = TRUE, col = "green")
curve_2D_fix_v(mylogit, v = 1, add = TRUE, col = "blue")
curve_2D_fix_v(mylogit, v = 1.2, add = TRUE, col = "cyan")
curve_2D_fix_v(mylogit, v = 0.4, add = TRUE, col = "yellow")
这是曲线:
讨论
在两个图中,我们看到 gender
(预测概率)和 v
(速度)之间的关系不是很强。在二维图中,几乎所有 v
的值都会产生相同的曲线。另一方面,press
(压力)是一个强大的影响。
回到你的模型:
> summary(mylogit)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.08326 4.45463 1.815 0.0696 .
press -0.02575 0.01618 -1.591 0.1115
v -0.15385 4.83824 -0.032 0.9746
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
可以看出v
一点都不显着!而严格来说,press
在 0.1 水平上也不显着。 所以这是一个非常弱的模型。我建议您删除变量 v
并再次创建模型,使用 press
作为唯一变量。
我想编写代码来绘制我的逻辑回归模型,即 "S" 形状的逻辑曲线。我有 两个独立的协变量 ,请问如何做到这一点?我附上了我的数据集和我的模型的代码。先感谢您。
239 0.72 1
324.6 0.83 1
331.8 0.95 1
334.3 0.83 1
259.7 0.89 1
212.3 0.88 1
204.7 0.65 1
253.86 0.75 1
258.94 0.85 1
329.66 0.95 0
469.68 1.46 0
459.74 1.11 0
293.2 0.64 0
297.88 0.98 0
267.9 0.82 0
374.1 1.29 0
333.62 0.74 0
dat <- read.table("data.txt")
colnames(dat)<-c("press","v","gender")
# logostic regression
dat$gender <- factor(dat$gender)
mylogit<- glm(gender~press+v,data=dat,family="binomial")
summary(mylogit)
######## the code below are irrelevant to making plot, ignore if you want
mylogit$fitted.values
newdat <- data.frame(t(c(300,0.1)))
colnames(newdat)<-c("press","v")
# this is your new dataset, we name it as "newdat"
pred <- predict(mylogit,newdata = newdat,type="response")
pred # the probability of being in class 1 will stored in this object
pred <- predict(mylogit,newdata = dat,type="response")
pred # the probability of being in class 1 will stored in this object
# accuracy
dat$pred <- 0
factor(dat$pred)
dat$pred[which(pred>0.5)] <- 1
table(dat$gender,dat$pred)
您有 2 个连续的非分类变量,因此逻辑曲线将是 3D 曲线。我给大家提供两种展示方式。
- 使用
persp
函数生成真实的3D平滑曲线; - 固定
v
为多个值,然后生成多个二维logistic曲线(你称之为"S"-shape曲线)。
3D曲线
press_grid <- seq(200, 480, by = 5)
v_grid <- seq(0.6, 1.5, by = 0.1)
newdat <- data.frame(press = rep(press_grid, times = length(v_grid)), v = rep(v_grid, each = length(press_grid)))
pred <- predict.glm(mylogit, newdata = newdat, type="response")
z <- matrix(pred, length(press_grid))
persp(press_grid, v_grid, z, xlab = "pressure", ylab = "velocity", zlab = "predicted probability", main = "logistic curve (3D)", theta = 30, phi = 20)
您需要先生成一个二维网格。 newdat
保存此网格,您可以执行 plot(newdat)
来查看此网格。然后通过调用 predict.glm(..., type = "response")
在此网格上进行预测。结果 pred
是一个向量。要绘制它,将其转换为矩阵 z
,然后调用 persp
进行 3D 绘图。 xlab
、ylab
、zlab
是三轴的标签。参数 theta
和 phi
用于调整视角。
在上面,press
和v
的边际网格是基于你的原始数据的范围:range(dat$press)
和range(dat$v)
。我们不做超出这个范围的预测。但即使在这个范围内,您也只有 17 个观测值。所以你仍然需要对情节持怀疑态度。
这是曲线:
二维曲线
这个玩具函数对于制作 2D 曲线很有用,v
固定为某个级别:
curve_2D_fix_v <- function(model, v = 1, press_grid = seq(200, 480, by = 5), add = FALSE, col = "black") {
newdat <- data.frame(press = press_grid, v = v)
pred <- predict.glm(model, newdat, type = "response")
if (add) lines(press_grid, pred, col = col) else {
plot(press_grid, pred, xlab = "pressure", ylab = "predicted probability", type = "l", col = col, main = "logistic curve (2D)")
abline(h = c(0, 0.5, 1), lty = 2, col = col)
}
}
如果add = FALSE
,则开启新的绘图window;虽然它是 TRUE
,但它绘制在之前的 window 上(但你有责任确保有这样一个 window!)二维图提供了更多信息,因为你可以添加一个0、0.5 和 1 处的水平线。
让我们开始吧:
curve_2D_fix_v(mylogit, v = 0.4, add = FALSE, col = "black")
curve_2D_fix_v(mylogit, v = 0.6, add = TRUE, col = "red")
curve_2D_fix_v(mylogit, v = 0.8, add = TRUE, col = "green")
curve_2D_fix_v(mylogit, v = 1, add = TRUE, col = "blue")
curve_2D_fix_v(mylogit, v = 1.2, add = TRUE, col = "cyan")
curve_2D_fix_v(mylogit, v = 0.4, add = TRUE, col = "yellow")
这是曲线:
讨论
在两个图中,我们看到 gender
(预测概率)和 v
(速度)之间的关系不是很强。在二维图中,几乎所有 v
的值都会产生相同的曲线。另一方面,press
(压力)是一个强大的影响。
回到你的模型:
> summary(mylogit)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.08326 4.45463 1.815 0.0696 .
press -0.02575 0.01618 -1.591 0.1115
v -0.15385 4.83824 -0.032 0.9746
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
可以看出v
一点都不显着!而严格来说,press
在 0.1 水平上也不显着。 所以这是一个非常弱的模型。我建议您删除变量 v
并再次创建模型,使用 press
作为唯一变量。