4
$\begingroup$

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?

$\endgroup$
2
  • 1
    $\begingroup$ 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. Try the manual calculation again but use more precise coefficients values. $\endgroup$ Dec 19, 2019 at 14:37
  • $\begingroup$ @COOLSerdash, wow, you are right! it was just because of not enough precise coefficients applied! now I will get my y as logit2prob(-1.270747) and this corresponds to expected value 0.2191294. Do you want to post this as an answer? $\endgroup$
    – maycca
    Dec 23, 2019 at 13:40

1 Answer 1

4
$\begingroup$

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$.

$\endgroup$
9
  • $\begingroup$ please, can you also specify how to get a coefficient for a reference categorical variable? i.e. for rank1? $\endgroup$
    – maycca
    Dec 30, 2019 at 7:41
  • $\begingroup$ @maycca I've added a short section. In essence, the coefficient of rank1 is the intercept. $\endgroup$ Dec 30, 2019 at 8:42
  • $\begingroup$ maybe I can just even skip the reference category from the prediction, as it always will be 0? $\endgroup$
    – maycca
    Dec 30, 2019 at 10:11
  • 1
    $\begingroup$ @maycca Maybe it's better to say that rank1 has no coefficient. That's probably less confusing. The coefficient of the ranks are the differences between the reference category, which is rank1 here and the other ranks. As the ranks are dummy-coded, rank1 is when all the dummy-variables for all other ranks are $0$. $\endgroup$ Jan 7, 2020 at 16:15
  • 1
    $\begingroup$ @maycca Yes, that's correct. $\endgroup$ Jan 7, 2020 at 16:40

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.