I have a data frame D1 in R with a dependent binary variable Response (0/1) and a set of covariates like age and gender. I want to know how "typical" everyone's Response is given their age and gender. I also have another, large data frame D2 with the same variables, which I can use to predict what is "typical". I fit a logistic regression model in D2, predicting Response from age and gender.

Now I want to take this model and use it to predict Response values for everyone in dataset D1, then compare these predicted values with the actual observed values, and finally get the residuals expressing this difference for every observation. In linear regression, this can be done simply by using the predict() function to calculate new values and then subtracting those from the observed values (like here). How can I do this in logistic regression, using a glm model in R?

The ultimate goal is to have residuals which I could then use in further analyses. These residuals would express how far away everyone is from what we would expect for their age and gender. I am aware there is a residuals() function available but that only seems to work for the data the original model was fit on.

See example below. I think I might be on the right track but I am not sure if my solution is doing what I need it to do.

# Create data frames D1 and D2
Response <- factor(rep(c(0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1), times = 4))
age <- rep(c(24, 31, 33, 89, 12, 67, 25, 91, 30, 56, 29, 80), times = 4)
gender <- factor(rep(c("M", "F"), times = 4, each = 6))
D1 <- data.frame(Response = Response[sample(48, 10)],
                 age = age[sample(48, 10)],
                 gender = gender[sample(48, 10)])
D2 <- data.frame(Response, age, gender)

# Create a model using D2
model <- glm(Response ~ gender + scale(age), family = binomial, data = D2, na.action = na.exclude)

# Predict values in D1 & extract residuals
observed <- as.numeric(D1$Response) - 1 # subtracting 1 from the value because as.numeric() turns a factor variable into 1/2 rather than 0/1
predicted <- predict(model, newdata = D1, type = "response") # Not sure if this is the type I want
residuals <- observed - predicted # Does this do what I want it to do?
glm_sanity_set <- cbind(D1, observed, predicted, residuals)

1 Answer 1


I haven't worked through the details of your code, but in outline it's a reasonable way to accomplish what you want. Quoting from the help page for predict.glm in R:

... for a default binomial model ... type = "response" gives the predicted probabilities.

So taking the difference between the observed 0/1 value and the predicted probability based on the new covariate values can be considered an estimate of how far from "typical" an individual's response is. The mean-square of all the differences is a form of the Brier score frequently used to evaluate the quality of a binary-outcome model, similar to how the mean-square error is used in ordinary least squares.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.