Get equation from glm coefficients: calculate y manually?

Question

I am trying to understand the math behind the glm(). Specifically, how to apply equation based on model predictors to calculate my y.pred? To fully understand it, I want to calculate my y manually, and apply this equation on the columns of dummy-coded predictors.

Let's have an example:

mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")

# data source: https://stats.idre.ucla.edu/r/dae/logit-regression/
# treat rank as factor   
mydata$rank <- factor(mydata$rank)

Define logistic regression as logit, using binomial family:

mylogit <- glm(admit ~ gre *gpa + rank, data = mydata, family = binomial("logit")) 

Ok, this is how my coefficients looks like:

Coefficients:
(Intercept)          gre          gpa        rank2        rank3        rank4  
 -13.608810     0.018344     3.652170    -0.721697    -1.343466    -1.606298  
    gre:gpa  
  -0.004719 

And this is the head of my input data:

  admit gre  gpa rank predict.admit
1     0 380 3.61    3     0.2191482
2     1 660 3.67    3     0.2948499
3     1 800 4.00    1     0.6398017
4     1 640 3.19    4     0.1886730
5     0 520 2.93    4     0.1029292
6     1 760 3.00    2     0.4522299

Following the logic that the predicted value is the sum of intercept + predictorsXbeta coefficients + beta*x1*x2 to include interactions, I get something like this:

 y = a + b1x1 + b2x2 + ... + bnxn + bx1*x2

  admit gre  gpa rank predict.admit
1     0 380 3.61    3     0.2191482

Combine the coefficients with the first row of my table:

  admit gre  gpa rank
1     0 380 3.61    3

 y = -13.608 +380*0.018 + 3.61*3.652 -1.343*1 - 0.004* 380*3.61

my y is -0.41448.

I thought that this is a logit(y) and as such I want to transform it back to probability value: function here

logit2prob <- function(logit){
  odds <- exp(logit)
  prob <- odds / (1 + odds)
  return(prob)
}

Let's convert my y.pred (-0.41448) to probability:

logit2prob(-0.41448)
[1] 0.3978384

But, if I run the predict.glm on my data, I get a different value:

mydata$predict.admit<-predict.glm(mylogit, mydata, type = "response")

admit gre  gpa rank predict.admit
1     0 380 3.61    3     0.2191482   # !! 0.2191482 is not equal to my calculation "-0.41448"??

So, I have in fact two questions:

1. How to calculate my y.predict manually, from estimated coefficients?
2. How to get y.predict for the rank = 1, i.e. the reference category, which is missing from the predicted values? is it simply 0?

My ultimate goal is to create a *vector of coefficients * that I can later on multiply by columns (binary for categorical values, qualitative for others) to predict my y manually. This requires convert the categorical values to binary dummy variables.

EDIT:

This is a full example to manually calculate my y using the dummy variables, where I skip the reference categories.

# get vector of coefficients
myCoef<- mylogit$coefficients

# get binary variables
mydata.bin <- fastDummies::dummy_cols(mydata,
                                      remove_first_dummy = TRUE)  # remove 

      admit gre  gpa rank rank_2 rank_3 rank_4
1     0 380 3.61    3      0      1      0
2     1 660 3.67    3      0      1      0
3     1 800 4.00    1      0      0      0
4     1 640 3.19    4      0      0      1
5     0 520 2.93    4      0      0      1
6     1 760 3.00    2      1      0      0

# remove not needed columns and add column for interaction
# I do this step to multiply vector of coefficients with individual columns, columnswise
md.bin<-subset(mydata.bin, select = -c(admit, rank))

# add new columns to include intercept and interactions between `gre` and `gpa`
md.bin$interc <- 1
md.bin$gre_gpa = md.bin$gre * md.bin$gpa

# order the dataframe as coefficients order
md.bin<- md.bin[c("interc", 
                  "gre",
                  "gpa",
                  "rank_2",
                  "rank_3",
                  "rank_4",
                  "gre_gpa" )]

# Multiply the vector of coefficients with dataframe
part.md.bin <- sweep(md.bin, 2, myCoef, "*")

# sum by rows and add intercept value
mydata.bin$pred.manual <- logit2prob(rowSums(part.md.bin))
mydata.bin$fitted <- mylogit$fitted.values 

It seems that my manually predicted values (pred.manual), and fitted values are the same (at some precision). I wonder if this means that my manual calculation is correct?

Answer 1

Your manual calculation is correct in terms of the coefficients, but you lose precision by either rounding incorrectly (i.e. the coefficient of the interaction) or by just not using enough numbers after the decimal point. Here is the manual calculation using 5 significant digits: $$ \mathrm{logit}(y) = -13.609 + 0.018344\cdot 380 + 3.6522\cdot 3.61 - 1.3435\cdot 1 -0.004719\cdot 380\cdot 3.61 = -1.270862 $$ The corresponding probability is: $$ \exp(-1.270862)/(1 + \exp(-1.270862)) = 0.2191 $$

Which is in good agreement with the more precise result from predict which uses the full precision available.

Edit

What if instead of the prediction for rank3 we'd like to predict for an individual with rank1 which is the reference value of the categorical variable rank? In the current dummy coding scheme, the coefficients for rank2-rank4 are the differences in the log-odds between the reference category and the respective level. Thus, the coefficient for the reference category is the intercept. Hence, the predicted log-odds in this situation are simply:

$$ \mathrm{logit}(y) = -13.609 + 0.018344\cdot 380 + 3.6522\cdot 3.61 -0.004719\cdot 380\cdot 3.61 = 0.0726378 $$ and the corresponding probability is $0.5181515$.