Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
#https://stats.stackexchange.com/questions/29044/plotting-confidence-intervals-for-the-predicted-probabilities-from-a-logistic-re
library(tidyverse)
library(magrittr)
set.seed(1234)
# create fake data on gambling. Does prob win depend on bid size?
mydat <- data.frame(
won=as.factor(sample(c(0, 1), 250, replace=TRUE)),
bid=runif(250, min=0, max=1000)
)
# logistic regression model:
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))
# new predictor values to use for prediction:
plotdat <- data.frame(bid=(0:1000))
# df with predictions, lower and upper limits of CIs:
preddat <- predict(mod1,
type = "link",
newdata=plotdat,
se.fit=TRUE) %>%
as.data.frame() %>%
mutate(bid = (0:1000),
# model object mod1 has a component called linkinv that
# is a function that inverts the link function of the GLM:
lower = mod1$family$linkinv(fit - 1.96*se.fit),
point.estimate = mod1$family$linkinv(fit),
upper = mod1$family$linkinv(fit + 1.96*se.fit))
# plotting with ggplot:
preddat %>% ggplot(aes(x = bid,
y = point.estimate)) +
geom_line(colour = "blue") +
geom_ribbon(aes(ymin = lower,
ymax = upper),
alpha = 0.5) +
scale_y_continuous(limits = c(0,1))
###https://stats.stackexchange.com/questions/167324/variance-covariance-matrix-of-logit-with-matrix-computation
library(Matrix)
library(sandwich)
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
mylogit <- glm(admit ~ gre + gpa, data = mydata, family = "binomial")
X <- as.matrix(cbind(1, mydata[,c('gre','gpa')]))
n <- nrow(X)
pi<-mylogit$fit
w<-pi*(1-pi)
v<-Diagonal(n, x = w)
var_b<-solve(t(X)%*%v%*%X)
var_b
x 3 Matrix of class "dgeMatrix"
[,1] [,2] [,3]
[1,] 1.1558251135 -2.818944e-04 -0.2825632388
[2,] -0.0002818944 1.118288e-06 -0.0001144821
[3,] -0.2825632388 -1.144821e-04 0.1021349767
vcov(mylogit)
(Intercept) gre gpa
(Intercept) 1.1558247051 -2.818942e-04 -0.2825631552
gre -0.0002818942 1.118287e-06 -0.0001144821
gpa -0.2825631552 -1.144821e-04 0.1021349526
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.