To calculate the "se.fit" and "resdiual.scale" that are the outputs from lm() function in R using python code

Question

i have a code snippet in R with simple lm() function with one dependent and one independent variable which is as follows.

X = ([149876.9876, 157853.421, 147822.3803, 147904.6639, 152625.6781, 147229.8083, 181202.081, 164499.6566, 171461.6586, 164309.3919])
Y = ([26212109.07, 28376408.76, 30559566.77, 26765176.65, 28206749.66, 27560521.33, 32713878.83, 31263763.7, 30812063.54, 30225631.6])

lmfit <- lm(formula = Data_df$Y ~ Data_df$X, data=Data_df)
lmpred <- predict(lmfit, newdata=Data_df, se.fit=TRUE, interval = "prediction")
print(lmpred) #prints out fit, se.fit, df, residual.scale

The output of the above code are 3 vectors 1.) fit 2.) se.fit 3.) df 4.) residual.scale

Please help me find the way to calculate se.fit and residual.scale in python. I m using statsmodels.ols to do the linear regression model. Below is the python code that i m using to build the linear regression.

import pandas as pd
import statsmodels.formula.api as smf
import numpy as np
ols_result = smf.ols(formula='Y ~ X', data=DATA_X_Y_OLS).fit()
ols_result.predict(data_x_values)

R output

$fit
        fit      lwr      upr
1  27594475 23262089 31926862
2  28768803 24486082 33051524
3  27291987 22943619 31640354
4  27304101 22956398 31651804
5  27999150 23686118 32312183
6  27204745 22851531 31557960
7  32206302 27951767 36460836
8  29747293 25490577 34004009
9  30772271 26527501 35017042
10 29719281 25462018 33976544


$se.fit
        1         2         3         4         5         6         7         8         9        10        
 578003.4  483363.7  605520.6  604399.0  542961.1  613642.7  420890.0  426036.9  397072.7  427318.3  

$df
[1] 24

$residual.scale
[1] 2017981

@akrun have given a sample of how the R output is. I should reproduce the same in python. — Anu Abraham, Feb 05 '19 at 17:37
Probably this can help: https://stackoverflow.com/q/38109501/4891738 — Zheyuan Li, Feb 05 '19 at 18:22

score 0 · Accepted Answer · edited Jul 23 '20 at 15:56

To find the fit, se.fit, df, residual.scale that are outputs of lm() function in R. Below is the python code to calculate the above mentioned 4 values

import statsmodels.formula.api as smf
import numpy as np
ols_result = smf.ols(formula='Y ~ X', data=DATA).fit()
fit = ols_result.predict(X_new) //predicted values ie, fit from lm()

covariance_matrix= ols_result.cov_params()
x = DATA['x'].values
xO = pd.DataFrame({"Constant":np.ones(len(x))}).join(pd.DataFrame(x)).values
x1 = np.dot(xO, COVARIANCE_MATRIX)
se_fit = np.sqrt(np.sum(x1 * xO,axis = 1)) //Standard error of the fitted values ie, se.fit in lm()
df = ols_result.df_resid //Degree of freedom ie, df in lm()

residual_scale = round(np.sqrt(np.dot(np.transpose(x), x)/df)) //Residual SD ie, Residual standard deviation

To calculate the "se.fit" and "resdiual.scale" that are the outputs from lm() function in R using python code

1 Answers1