# Difference Between GLM predict Output in R vs. GLM predict Output in Python (statsmodels) Logistic Regression

I am just working on trying to port some R script to Python. I am fairly new to the Python language. I have been going through the R script and converting it. But there is a part that is stumping me, where I ran a logistic regression in statsmodels GLM, but when I use the predict method afterwards, the results I am getting from statsmodels is not comparable to what I get running R.

I just need some help in figuring out how Python is running, to lead to the differences. To make it easier for reproducibility, I will include both my R and Python scripts. I will post my R scripts first, followed by the Python along with the data.

# R Script

# the dataframe
dfR = data.frame('DV' = c(1, 1, 0, 0, 1, 0, 1, 0, 0, 1),
'IV1' = c(3.58517, 2.07391, 0.55865, 1.08131, 0.15004, 2.9959, 1.29169, 1.50527, 0.81436, -1.90697),
'IV2' = c(3.73935, 0.53802, 1.73509, 3.01428, 2.7058, 1.97086, 3.98669, 2.67305, 3.18234, 1.39142),
'IV3' = c(0.8255, 4.8461, 4.3388, 2.565, 5.6754, 5.4174, 7.666, -0.1009, 0.4374, 3.6268)
)
print(dfR)

Y <- dfR[,1]
X <- dfR[,2:length(dfR[1,])]
XS <- scale(X)

XS.svd<-svd(XS)
Q<-XS.svd$$v P<-XS.svd$$u
Z<-P%*%t(Q)
ZS <- scale(Z)

Lambda<-solve(t(ZS)%*%ZS)%*%t(ZS)%*%XS

logrfit <- glm(unlist(Y) ~ ZS, family=binomial)
summary(logrfit)

unstCoefs<-coef(logrfit)
unstCoefs
b<-unstCoefs[2:length(unstCoefs)]
b

# **where there is a difference between R and Python**
LpredY<-predict(logrfit, newdata=dfR,type="response")
print(LpredY)


Below is the Python code, it follows the R code exactly, with just minor variations to be consistent with how Python outputs stuff

# Python Code

import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from numpy import linalg
import statsmodels.api as sm
import statsmodels.formula.api as smf

# the dataframe
dfP = pd.DataFrame(
{'DV': [1, 1, 0, 0, 1, 0, 1, 0, 0, 1],
'IV1': [3.58517, 2.07391, 0.55865, 1.08131, 0.15004, 2.9959, 1.29169, 1.50527,
0.81436, -1.90697],
'IV2': [3.73935, 0.53802, 1.73509, 3.01428, 2.7058, 1.97086, 3.98669, 2.67305,
3.18234, 1.39142],
'IV3': [0.8255, 4.8461, 4.3388, 2.565, 5.6754, 5.4174, 7.666, -0.1009, 0.4374, 3.6268]
}
)
print(dfP)

Y = dfP.iloc[:, 0]
X = dfP.iloc[:, 1:]
XS = StandardScaler().fit_transform(X)

P, _, Q = np.linalg.svd(XS, full_matrices=False)
Z = P @ Q
ZS = StandardScaler().fit_transform(Z)

Lambda = linalg.inv(ZS.T @ ZS) @ ZS.T @ XS

dfZS = pd.DataFrame(ZS, columns=X.columns)
dfXY = pd.concat([Y, dfZS], axis=1)

ZS_columns = list(dfXY.iloc[:, 1:].columns)
ZS2_columns = ' + '.join(str(i) for i in ZS_columns)

formula = 'Y ~ {}'.format(ZS2_columns)
mod = smf.glm(formula=formula, data=dfXY, family=sm.families.Binomial()).fit()
print(mod.summary())

b = mod.params[1:]

# **where there is a difference between R and Python**
LpredY = mod.predict(dfP)
print(LpredY)


The only part that is not consistent with the output from R is the LpredY variable. Just need some assistance in figuring out the discrepancies.

• Hi, statsmodels doesn't include intercept terms by default; try adding a column of ones to your data. Commented Aug 11, 2022 at 15:24
• I used the formula approach in statsmodels, which I think, defaults to inclusion of the intercept
– GSA
Commented Aug 11, 2022 at 15:37
• Oh TIL (as long as your output has a row that says "intercept" :)). Commented Aug 11, 2022 at 15:48

Two issues:

1. (minor): Python's StandardScaler differs slightly from R's scale. While scale in R uses $$1 / (n-1)$$ when calculating the standard deviation (the unbiased estimator), StandardScaler uses $$1/n$$ (max. likelihood). That's easy to fix: In Python, after defining XS and ZS, append XS /= np.sqrt(10/9) and ZS /= np.sqrt(10/9), respectively.

2. (major): In your R code, you fit the model using a vector (Y) and a matrix (ZS), but you run predict on a data frame (dfR). I'm not familiar with the implementation details in R, but I know that this is not the way to do it. The correct way would be something like:

    dfZS <- data.frame(Y, ZS)
colnames(dfZS) <- c('DV', 'IV1', 'IV2', 'IV3')
logrfit <- glm(DV ~ IV1 + IV2 + IV3, data=dfZS, family=binomial)


The output in R is then:

> print(LpredY)
1         2         3         4         5         6         7         8         9
0.6205818 0.9898041 0.9897777 0.9420639 0.9979161 0.9941012 0.9996996 0.4752596 0.6651955
10
0.9877343


and in Python:

In [182]: print(LpredY)
0    0.620582
1    0.989804
2    0.989778
3    0.942064
4    0.997916
5    0.994101
6    0.999700
7    0.475260
8    0.665195
9    0.987734
dtype: float64


After much reading, I came upon the solution. It turns out that the following script in R GLM, family=binomial

predict(logrfit, newdata=dfR, type="response")


has its counterpart in Python statsmodel GLM, family=sm.families.Binomial() as

model.mu


This means that in the question I posted earlier for this R script:

predict(logrfit, newdata=dfR, type="response")


The equivalent Python form is not mod.predict(dfP), but rather,

mod.mu


where mod is the name I gave my GLM object

Below are some links to the information on that: yhat = res.mu and mu: ndarray

• predict used in-sample and mu should be the same numbers. Maybe dfR, dfP is not the appropriately transformed design matrix for prediction. Commented Aug 12, 2022 at 18:19