CI for predicted probabilities

[MarieAgosta]

Hi all,

So we manually did the predicted probabilities and now have to estimate the confidence intervals for our predicted probabilities… but the predict argument does not work for us.

We’ve searched around for ways to go about it and have yet to find anything that works. Do you have any suggestions? Below is the code, the bolded bit at the end is the new trouble spot.

Many thanks,

Marie

#manually collecting probabilities

####Black Millennial Woman
#1-probability for below-high-school educated black employed woman who didn't vote in 2008
prob_genY_ed1 <- (exp(-0.3448 - 0.3448 + 0.1249 + 0.6389))/(1+exp(-0.3448 - 0.3448 + 0.1249 + 0.6389)) 
#2-probability for high-school educated black employed woman who didn't vote in 2008
prob_genY_ed2 <- (exp(-0.3448 + 0.5263 + 0.1249 + 0.6389))/(1+exp(-0.3448 + 0.5263 + 0.1249 + 0.6389)) 
#3-probability for some post-high-school educated black employed woman who didn't vote in 2008
prob_genY_ed3 <- (exp(-0.3448 + 0.7955 + 0.1249 + 0.6389))/(1+exp(-0.3448 + 0.7955 + 0.1249 + 0.6389)) 
#4-probability for bachelor-educated black employed woman who didn't vote in 2008
prob_genY_ed4 <- (exp(-0.3448 + 0.8990 + 0.1249 + 0.6389))/(1+exp(-0.3448 + 0.8990 + 0.1249 + 0.6389)) 
#5-probability for graduate-educated black employed woman who didn't vote in 2008
prob_genY_ed5 <- (exp(-0.3448 + 1.5099 + 0.1249 + 0.6389))/(1+exp(-0.3448 + 1.5099 + 0.1249 + 0.6389)) 


####Black Baby Boomer Woman
#1-probability for below-high-school educated black employed woman who didn't vote in 2008
prob_boomer_ed1 <- (exp(-0.87481 - 0.87481 + 0.22553 -0.76681))/(1+exp(-0.87481 - 0.87481 + 0.22553 -0.7668)) 
#2-probability for high-school educated black employed woman who didn't vote in 2008
prob_boomer_ed2 <- (exp(-0.87481 + 0.81340 + 0.22553 -0.76681))/(1+exp(-0.87481 + 0.81340 + 0.22553 -0.76681)) 
#3-probability for some post-high-school educated black employed woman who didn't vote in 2008
prob_boomer_ed3 <- (exp(-0.87481 + 1.12412 + 0.22553 -0.76681))/(1+exp(-0.87481 + 1.12412 + 0.22553 -0.76681)) 
#4-probability for bachelor-educated black employed woman who didn't vote in 2008
prob_boomer_ed4 <- (exp(-0.87481 + 1.69306 + 0.22553 -0.76681))/(1+exp(-0.87481 + 1.69306 + 0.22553 -0.76681)) 
#5-probability for graduate-educated black employed woman who didn't vote in 2008
prob_boomer_ed5 <- (exp(-0.87481 + 1.54687 + 0.22553 -0.76681))/(1+exp(-0.87481 + 1.54687 + 0.22553 -0.76681)) 

#creating a dataframe of predicted variables for plotting and other uses
blackfem_genY_ed <- c(prob_genY_ed1, prob_genY_ed2, prob_genY_ed3, prob_genY_ed4, prob_genY_ed5)
blackfem_boomer_ed <- c(prob_boomer_ed1, prob_boomer_ed2, prob_boomer_ed3, prob_boomer_ed4, prob_boomer_ed5)
education <- cbind(blackfem_genY_ed, blackfem_boomer_ed)

#creating a data frame of fitted values for a hypothetical black woman

blackfem_genY <- with(anes_genY, data.frame(female = 1, race = 2, vote_2008 = 0, unemployed = 0, education = factor(1:5))) 

#this is an alternative we've played around with that includes the weighting columns
#used by the survey package as well
blackfem_genY <- with(anes_genY, data.frame(caseID = factor(1:5), weights = .506, psu = 1, strata = 6, female = 1, race = 2, vote_2008 = 0, unemployed = 0, education = factor(1:5))) 

#predict probability point estimates for each fitted value.

blackfem_genY$predicted_prob <- prediction(M_genY, data = blackfem_genY, type = "response")

blackfem_genY$predicted_prob <- predict(M_genY, newdata = blackfem_genY, type = "response", se = TRUE)
#another alternative
blackfem_genY$predicted_prob <- predict(M_genY, newdata = blackfem_genY, type = "response", 
                                        design = ANESdesign_genY)

#Now to plot this with confidence intervals.
#First, add the confidence intervals

blackfem_genY <- cbind(blackfem_genY, predict(M_genY, newdata = blackfem_genY, type="link", se=TRUE))
blackfem_genY <- within(blackfem_genY, {
  prob <- plogis(fit)
  LL <- plogis(fit - (1.96 * se.fit))
  UL <- plogis(fit + (1.96 * se.fit))
})

[christophergandrud]

Hi @MarieAgosta

Some small notes to make it easier for others to help:

• Provide a link to the data so that the issue is reproducible.

• More clearly specify what you mean by it "doesn't work".

• It's good practice not to overload the question. Put only the minimal amount of code and data needed to reproduce the problem.

Re predicted probabilities: you are repeating a lot of code. Whenever you find yourself repeating a lot of code, it's time to create a function, e.g. something like:

# Note assumes a model with 3 explanatory variables
predicted_fit <- function(alpa, beta1, beta2, beta3, x1, x2, x3){
    fit <- alpha + beta1 * x1 + beta2 * x2 + beta3 * x3
    pred_prob <- (exp(fit)) / (1 + exp(fit)) 
    return(pred_prob)
}

You then enter the values like: sample_fitted <- predicted_fit(alpha = 1, beta1 = 0.434 and so on.

To really cut down on the code repetition you could then stick this in a loop, but don't worry about that now.

Confidence intervals

I'm a little confused about what is going on here. It seems like you are manually finding the predicted probabilities and then also using prediction to find them? Why do both?

For the confidence intervals, I'm confused about your use of plogis.

Wouldn't you instead want find the fitted values, then use these to come up with the confidence intervals? You could combine these steps with the above function, something like:

# Note assumes a model with 3 explanatory variables
predicted_fit <- function(alpha, beta1, beta2, beta3, x1, x2, x3, se){
    fit <- alpha + beta1 * x1 + beta2 * x2 + beta3 * x3
    lb <- fit - (1.96 * se)
    ub <- fit + (1.96 * se) 
    predictions <- c(lb, fit, ub)

    probs <- function(x) {(exp(x)) / (1 + exp(x) }
    predictions <- probs(predictions)
    return(predictions)
}

You could then rbind multiple results of this into a data frame and plot.

Note that I haven't actually run this code, but it probably is a good start.

[MarieAgosta]

Thanks professor and sorry, the predict argument is in the old code, we just haven't removed it yet. I shouldn't have included it.

[christophergandrud]

No worries. Let me know how it goes.

[MarieAgosta]

Thank you professor. We tried out the code, we adapted the code from what you wrote but we're struggling with the standard errors. There are five coefficients in our model so there will be five standard errors. How would we account for all of the standard errors?

Attempt using suggestions from prof

predicted_fit <- function(alpha, beta1, beta2, beta3, beta4, beta5, x1, x2, x3, x4, x5, se1, se2, se3, se4, se5)
{
fit <- alpha + beta1 * x1 + beta2 * x2 + beta3 * x3 + beta4 x4 + beta5 x5
lb <- fit - (1.96 * se)
ub <- fit + (1.96 * se)
predictions <- c(lb, fit, ub)

probs <- function(x) {(exp(x)) / (1 + exp(x)) }
predictions <- probs(predictions)
return(predictions) 

}

[christophergandrud]

Good point, sorry I neglected that earlier.

predicted_fit <- function(alpha, beta1, beta2, beta3, beta4, beta5, 
    x1, x2, x3, x4, x5, 
    se1, se2, se3, se4, se5)
{
lb_fun <- function(coef, se){
    lb <- coef - (1.96 * se)
    return(lb)
}

ub_fun <- function(coef, se){
    ub <- coef + (1.96 * se)
    return(ub)
}

fit <- alpha + beta1 * x1 + beta2 * x2 + beta3 * x3 + beta4 x4 + beta5 x5
lb <- alpha + lb_fun(beta1, se1) * x1  + lb_fun(beta2, se2) * x2  + lb_fun(beta3, se3) * x3
ub <- alpha + ub_fun(beta1, se1) * x1  + ub_fun(beta2, se2) * x2  + ub_fun(beta3, se3) * x3
predictions <- c(lb, fit, ub)

probs <- function(x) {(exp(x)) / (1 + exp(x)) }
predictions <- probs(predictions)
return(predictions) 
}

This code is really nasty and really the best way to do this would be to hack into the fitted model object and extract these programatically, but I don't have much more time today. Give this a shot and see what happens.

[MarieAgosta]

Thank you!!

[kardishcj]

Actually one quick follow-up question. For the upper and lower bounds, should the code stop at the 3rd beta or go all the way to the 5th (or however many coefficients there are in the model)? It does appear to work, though, yes.
fit <- alpha + beta1 * x1 + beta2 * x2 + beta3 * x3 + beta4 * x4 + beta5 * x5
lb <- alpha + lb_fun(beta1, se1) * x1 + lb_fun(beta2, se2) * x2 + lb_fun(beta3, se3) * x3
ub <- alpha + ub_fun(beta1, se1) * x1 + ub_fun(beta2, se2) * x2 + ub_fun(beta3, se3) * x3
predictions <- c(lb, fit, ub)

[christophergandrud]

Sorry, right it should include all of the coefficient point estimates beta I just did 3 for expediency.

[kardishcj]

Right, I figured. Is it possible that there's an error in the confidence-interval calculation? It's giving very wide confidence bands, which of course look odd when plotting. Or is that simply a function of having very high standard errors and imprecise estimates? Looking at the model output, they certainly don't seem that high at least.

[christophergandrud]

Can you given an example, including both the coefficient estimates/standard errors and predicted probabilities?

Remember that it is difficult to mentally relate these two as they are are different scales--hence the need to do predicted probabilities.

Also, stepping back, in class I did cover a simulation approach to estimating uncertainty that is especially useful when it is numerically tricky to mathematically come up with the uncertainty bounds. Do take a look.

[kardishcj]

I can take a look at the simulation approach. Below is a text-version stargazer table of the model across the four generations, followed by predicted probabilities (w/ upper and lower bounds) for a black, female Millennial of moderate income. I apologize for how the table looks but it doesn't really copy-paste well. Parentheses are standard errors.

                  Dependent variable:         
         -------------------------------------
             Probability of voting in 2012    
          Gen Y    Gen X    Boomers   Silent  
           (1)            (2)         (3)            (4)   

---

education 1.317*** 1.133*** 1.377*** 1.404***
(0.134) (0.118) (0.120) (0.192)

female 1.112*** 1.360*** 1.206*** 0.484
(0.235) (0.247) (0.253) (0.481)

vote_2008 2.956*** 6.727*** 14.474*** 23.935***
(0.248) (0.271) (0.254) (0.625)

black 2.640*** 1.567*** 0.482 6.277***
(0.393) (0.426) (0.357) (0.869)

hispanic 1.878*** 0.680** 0.883** 0.262
(0.352) (0.298) (0.390) (0.832)

income 1.217*** 1.161*** 1.196*** 1.665***
(0.079) (0.087) (0.069) (0.215)

Constant 0.240 0.378 0.224 0.144
(0.413) (0.495) (0.444) (0.812)

Predicted probabilities for black, female Millennial at different levels of education and moderate income

Below high school
lower: 0.1458708
predicted: 0.4814261
upper: 0.4199645

High school
lower: 0.1476030
predicted: 0.5501755
upper: 0.5534615

Some post-high
lower: 0.1493521
predicted: 0.6170605
upper: 0.6796692

Bachelor
lower: 0.1511182
predicted: 0.6797885
upper: 0.7841197

Graduate
lower: 0.1529015
predicted: 0.7366269
upper: 0.8614545

[christophergandrud]

Those don't look wildly off. We might be missing something, but I would check by comparing them to the simulated intervals.

[kardishcj]

Strange, but my Zelig simulations look much different than my predicted probabilities (shown below as plotted on ggplot).

This is the ggplot of the predictions for a black female of average income at various levels of education (top line is Millennial, bottom is Baby Boomer, still need to work on the aesthetics, I know)

Then below is a Millennial followed by Baby Boomer simulation using the same characteristics as the ggplot above.

Millennial

Baby Boomer

These Zelig simulations look quite different from the ggplot above. And there's not much of a distinction between the two simulations, which is another stark contrast with the ggplot.

[christophergandrud]

What are the solid lines in the ggplot one?

[christophergandrud]

btw, I would go with the Zelig predictions (we may have messed up the formula somewhere along the way) and they are much more straightforward to interpret.