代码之家 › 专栏 › 技术社区 › Dyllan

一阶差系数图

probability plot r

Dyllan · 技术社区 · 11 年前

我在r中运行一个逻辑模型,我试图用系数图表示自变量的概率差异。具体来说,我想通过将感兴趣的变量从其最小值转移到其最大值(同时保持其他变量的均值或模式)来创建概率差异。

在所附的图片中,我希望我的图表看起来与上半部分相似。 enter image description here

我运行了以下代码:

 mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
 mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family =
 "binomial")

然后,我计算了每个变量的最小值和最大值的预测概率,并将两者相减。我在区间的上下限重复了这个过程。附件是我的代码

plotdat <- data.frame(gre=c(.220, 800), gpa=mean(mydata$gpa, na.rm=TRUE), rank=c(2) ) 
preddat <- predict(mylogit, newdata=plotdat, se.fit=TRUE)

Grebeta<-(exp(preddat$fit[2])/(1+exp(preddat$fit[2])))-(exp(preddat$fit[1])/(1+exp(preddat$fit[1])))
Gremin<-(exp(preddat$fit[2]+1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]+1.96*preddat$se.fit[2])))-exp(preddat$fit[1]+1.96*preddat$se.fit[1])/(1+exp(preddat$fit[1]+1.96*preddat$se.fit[1]))
Gremax<-exp(preddat$fit[2]-1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]-1.96*preddat$se.fit[2]))-exp(preddat$fit[1]-1.96*preddat$se.fit[1])/(1+exp(preddat$fit[1]-1.96*preddat$se.fit[1]))


plotdat <- data.frame(gpa=c(2.26, 4), gre=mean(mydata$gre, na.rm=TRUE), rank=c(2) )
preddat <- predict(mylogit, newdata=plotdat, se.fit=TRUE)

GPAbeta<-(exp(preddat$fit[2])/(1+exp(preddat$fit[2])))-(exp(preddat$fit[1])/(1+exp(preddat$fit[1])))
GPAmin<-(exp(preddat$fit[2]+1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]+1.96*preddat$se.fit[2])))-exp(preddat$fit[1]+1.96*preddat$se.fit[1])/(1+exp(preddat$fit[1]+1.96*preddat$se.fit[1]))
GPAmax<-exp(preddat$fit[2]-1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]-1.96*preddat$se.fit[2]))-exp(preddat$fit[1]-1.96*preddat$se.fit[1])/(1+exp(preddat$fit[1]-1.96*preddat$se.fit[1]))


plotdat <- data.frame(rank=c(4, 1), gre=mean(mydata$gre, na.rm=TRUE), gpa=mean(mydata$gpa, na.rm=TRUE ))
preddat <- predict(mylogit, newdata=plotdat, se.fit=TRUE)

Rankbeta<-(exp(preddat$fit[2])/(1+exp(preddat$fit[2])))-(exp(preddat$fit[1])/(1+exp(preddat$fit[1])))
Rankmin<-(exp(preddat$fit[2]+1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]+1.96*preddat$se.fit[2])))-exp(preddat$fit[1]+1.96*preddat$se.fit[1])/(1+exp(preddat$fit[1]+1.96*preddat$se.fit[1]))
Rankmax<-exp(preddat$fit[2]-1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]-1.96*preddat$se.fit[2]))-exp(preddat$fit[1]-1.96*preddat$se.fit[1])/(1+exp(preddat$fit[1]-1.96*preddat$se.fit[1]))

然后,我创建了三个包含概率和频带差异的向量。附件是我的代码:

se.max<- c(Gremax   , GPAmax  , Rankmax  )
coef.vec<- c( Grebeta  ,GPAbeta  , Rankbeta ) 
se.min<-c(Gremin , GPAmin, Rankmin)



var.names <- c("gre", "gpa", "rank")

最后,我绘制了图表。

y.axis <- c(length(coef.vec):1)

par(mar=c(2, 13, 0, 0))

plot(coef.vec, y.axis, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.2,  xlim = c(-2,2.5), xaxs = "r", main = "")

segments(se.max, y.axis,se.min, y.axis, lwd =  1.5)

axis(1, at = seq(-1,1,by=.25), labels = NA, tick = T,cex.axis = 1.2, mgp = c(2,.7,0))
axis(1, at = seq(-1,1,by=.5), labels =  c(-1,  -.5,  0, .5,1), tick = T,cex.axis = 1.2, mgp = c(2,.7,0))

axis(2, at = y.axis, label = var.names, las = 1, tick = T, ,mgp = c(2,.6,0), cex.axis = 1.2)
segments(0,0,0,17,lty=2)

然而,我无法绘制置信区间。以下是我的最终输出。

enter image description here

看来我的信心指数不会有变化。如果有人能提供帮助并指出我的计算或代码中的错误,我将非常感激。

1 回复 | 直到 11 年前

Dyllan 11 年前

plotdat <- data.frame(gre=c(.220, 800), gpa=mean(mydata$gpa, na.rm=TRUE), rank=c(2) ) 
preddat <- predict(mylogit, newdata=plotdat, se.fit=TRUE)

#GRE High
GREbetahigh<-(exp(preddat$fit[2])/(1+exp(preddat$fit[2])))
GREminhigh<-(exp(preddat$fit[2]+1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]+1.96*preddat$se.fit[2])))
GREmaxhigh<-exp(preddat$fit[2]-1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]-1.96*preddat$se.fit[2]))


#GRE low
GREbetalow<-(exp(preddat$fit[1])/(1+exp(preddat$fit[1])))
GREminlow<-(exp(preddat$fit[1]+1.96*preddat$se.fit[2])/(1+exp(preddat$fit[1]+1.96*preddat$se.fit[1])))
GREmaxlow<-exp(preddat$fit[1]-1.96*preddat$se.fit[2])/(1+exp(preddat$fit[1]-1.96*preddat$se.fit[1]))

#GRE Diff
GREbeta.diff<-GREbetahigh-GREbetalow
GREmax.diff<-GREmaxhigh-GREmaxlow
GREmin.diff<-GREminhigh-GREminlow

#GPA
plotdat <- data.frame(gpa=c(2.26, 4), gre=mean(mydata$gre, na.rm=TRUE), rank=c(2) )
preddat <- predict(mylogit, newdata=plotdat, se.fit=TRUE)

#GPA high
GPAbetahigh<-(exp(preddat$fit[2])/(1+exp(preddat$fit[2])))
GPAminhigh<-(exp(preddat$fit[2]+1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]+1.96*preddat$se.fit[2])))
GPAmaxhigh<-exp(preddat$fit[2]-1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]-1.96*preddat$se.fit[2]))

#GPA low
GPAbetalow<-(exp(preddat$fit[1])/(1+exp(preddat$fit[1])))
GPAminlow<-(exp(preddat$fit[1]+1.96*preddat$se.fit[2])/(1+exp(preddat$fit[1]+1.96*preddat$se.fit[1])))
GPAmaxlow<-exp(preddat$fit[1]-1.96*preddat$se.fit[2])/(1+exp(preddat$fit[1]-1.96*preddat$se.fit[1]))

#GPA Diff
GPAbeta.diff<-GPAbetahigh-GPAbetalow
GPAmax.diff<-GPAmaxhigh-GPAmaxlow
GPAmin.diff<-GPAminhigh-GPAminlow

#Rank

plotdat <- data.frame(rank=c(4, 1), gre=mean(mydata$gre, na.rm=TRUE), gpa=mean(mydata$gpa, na.rm=TRUE ))
preddat <- predict(mylogit, newdata=plotdat, se.fit=TRUE)

#Rank high
Rankbetahigh<-(exp(preddat$fit[2])/(1+exp(preddat$fit[2])))
Rankminhigh<-(exp(preddat$fit[2]+1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]+1.96*preddat$se.fit[2])))
Rankmaxhigh<-exp(preddat$fit[2]-1.96*preddat$se.fit[2])/(1+exp(preddat$fit[2]-1.96*preddat$se.fit[2]))

#Rank Low
Rankbetalow<-(exp(preddat$fit[1])/(1+exp(preddat$fit[1])))
Rankminlow<-(exp(preddat$fit[1]+1.96*preddat$se.fit[1])/(1+exp(preddat$fit[1]+1.96*preddat$se.fit[1])))
Rankmaxlow<-exp(preddat$fit[1]-1.96*preddat$se.fit[1])/(1+exp(preddat$fit[1]-1.96*preddat$se.fit[1]))


#Rank Diff
Rankbeta.diff<-Rankbetahigh-Rankbetalow
Rankmax.diff<-Rankmaxhigh-Rankmaxlow
Rankmin.diff<-Rankminhigh-Rankminlow

#Graph
se.max<- c(GREmax.diff   , GPAmax.diff, Rankmax.diff)
coef.vec<- c( GREbeta.diff , GPAbeta.diff, Rankbeta.diff)
se.min<-c(GREmin.diff , GPAmin.diff, Rankmin.diff)

var.names <- c("gre", "gpa", "rank")

y.axis <- c(length(coef.vec):1)

par(mar=c(2, 13, 0, 0))


plot(y.axis, coef.vec, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.2,  ylim = c(-1,1), xlim=c(1,3.3), xaxs = "r", main = "")
segments(y.axis, se.max,y.axis, se.min, lwd =  1.5)

axis(2, at = seq(-1,1,by=.25), labels = NA, tick = T,cex.axis = 1.2, mgp = c(2,.7,0))
axis(2, at = seq(-1,1,by=.5), labels =  c(-1,  -.5,  0, .5, 1), tick = T,cex.axis = 1.2, mgp = c(2,.7,0))

axis(1, at = y.axis, label = var.names, las = 1, tick = T, ,mgp = c(2,.6,0), cex.axis = 1.2)
segments(1,0,3.3,0,lty=2)