Skip to content

Instantly share code, notes, and snippets.

# adamwespiser/conf intervals for point-wise prediction Last active Oct 18, 2018

 #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
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.