7.1 Using the response functions to generate quantities of interest

Recall, in linear regression, to get our estimated values \(\hat Y\) we said \(\hat Y = X\hat\beta\).

  • In glm’s, we can do the same to get our estimated values on the scale of the linear predictor \(\hat \eta = X\hat\beta\).
  • We then use our \(Link^{-1}\) response function to transform these values into the quantity of interest.
    • E.g., in logistic regression we want \(\hat{\pi} = \frac{\exp(X\hat\beta)}{1 + \exp(X\hat\beta)}\).
    • E.g., in probit regression we want \(\hat{\pi} = \Phi(X \hat \beta)\).
    • These represent the predicted probability of \(Y_i = 1\) given our coefficient estimates and designated values of the covariates

Let’s use a subset of the MIDs mids.txt data available here.

This dataset has variables related to whether a dyad of states is engaged in a militarized interstate dispute between the two countries in a given year. The variable that will be our outcome of interest is Conflict which takes the values 0 or 1. We will also look at the relationship between a few independent variables and the propensity for conflict. Data are at Dyad Level.

  • whether the pair of countries
    • include a major power (MajorPower, 1=yes, 0=otherwise),
    • are contiguous ( Contiguity, 1=yes, 0=otherwise),
    • are allies (Allies, 1=yes, 0=otherwise),
    • and/or have similar foreign policy portfolios (ForeignPolicy, 1=yes, 0=otherwise)
  • BalanceofPower: balance of military power
  • YearsSince: the number of years since the last dispute
## Load data
mids <- read.table("https://raw.githubusercontent.com/ktmccabe/teachingdata/main/midsshort.txt")

table(mids$Conflict)

    0     1 
99657   343

As can be seen from the table above, conflicts (fortunately) are relatively rare in our data. This means are predicted probabilities are likely going to be pretty small in this example.

We will run a logistic regression with a few covariates.

out.logit <-glm(Conflict ~  MajorPower + Contiguity + Allies + ForeignPolicy +
    BalanceOfPower + YearsSince, 
                family = binomial(link = "logit"), data = mids)

Our logistic regression equation is:

  • \(\log \frac{\pi_i}{1-\pi_i} = \mathbf{x_i'}\beta\), or alternatively
  • \(Pr(Y_i = 1 | \mathbf{x}_i) = logit^{-1}(\mathbf{x}_i'\beta) = \frac{exp(\mathbf{x_i'}\beta)}{1 + exp(\mathbf{x_i'}\beta)}\)

Our coefficients are in the “logit” aka log-odds scale of the linear predictor, so we use the response function to put them into probability estimates.

For logistic regression, we can generate predicted probabilities for each \(Y_i\) using the

  1. predict(model, type="response") function, using the
  2. plogis function with \(X\hat \beta\), or
    • Question: What would be the equivalent function for probit?
  3. manually writing down the response function \(\frac{exp(\mathbf{x_i'}\beta)}{1 + exp(\mathbf{x_i'}\beta)}\).

We use the predict function exactly the same way as before, but the key argument we need to specify is type. If you have type = link this generates answers that are still on the log-odds linear predictor scale. It is switching this to response that goes into probability in the logit case or probit case.

Example

First, let’s just generate predicted probabilities for each observation in our data, without specifiying any designated values for \(X\)– just keeping the values as they are in the data aka “as observed.”

## Method with predict()
## When you don't specify newdata, R assumes you want the X data from the model
pp <- predict(out.logit, type = "response")

## Manual way #1
X <- model.matrix(out.logit)
bh <- coef(out.logit)
pplogis <- plogis(X %*% bh)

## Manual way # 2
ppexp <- exp(X %*% bh)/(1 + exp(X %*% bh))

## Compare the first five rows of each to see they are the same
cbind(pplogis, ppexp, pp)[1:5,]
                                           pp
78627  0.0004204501 0.0004204501 0.0004204501
295818 0.0001035526 0.0001035526 0.0001035526
251841 0.0006211178 0.0006211178 0.0006211178
98068  0.0004270812 0.0004270812 0.0004270812
209797 0.0001005396 0.0001005396 0.0001005396

The code above generates a predicted probability associated with each observation in the data. This is similar to generating a fitted value \(\hat y\) for each observation in OLS.