1
$\begingroup$

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.

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

2 Answers 2

2
$\begingroup$

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
$\endgroup$
0
$\begingroup$

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

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

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.