Simple linear regression of two variables

Question

This is the code I built to implement a simple linear regression to check the impact of variable A on B. I am not interested in predicting that's why the statistical approach is important for me to focus on the descriptive analytics such as the statsmodels summary and comparing RMSE of training set and test set to see if the model did a good job or not (I'm not looking to compare Dependent variable and RMSE to check prediction quality for example). Is there anything fundamentally wrong with this approach to a simple Linear Regression?

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy import stats

# Import Excel File
data = pd.read_excel("C:\\Users\\AchourAh\\Desktop\\Simple_Linear_Regression\\SP Level Simple Linear Regression\\PL32_PMM_03_09_2018_SP_Level.xlsx",'Sheet1') #Import Excel file

# Replace null values of the whole dataset with 0
data1 = data.fillna(0)
print(data1)

# Extraction of the independent and dependent variable
X = data1.iloc[0:len(data1),1].values.reshape(-1, 1) #Extract the column of  the COPCOR SP we are going to check its impact
Y = data1.iloc[0:len(data1),2].values.reshape(-1, 1) #Extract the column of the PAUS SP

# Data Splitting to train and test set
from sklearn.model_selection import train_test_split
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size =0.25,random_state=42)


# Training the model and Evaluation of the Model
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
import math
from sklearn import model_selection
from sklearn.model_selection import KFold
lm = LinearRegression() #create an lm object of LinearRegression Class
lm.fit(X_train, Y_train)  #train our LinearRegression model using the   training set of data - dependent and independent variables as parameters.  Teaching lm that Y_train values are all corresponding to X_train.
kf = KFold(n_splits=6, random_state=None)
for train_index, test_index in kf.split(X_train):
print("Train:", train_index, "Validation:",test_index)
X_train1, X_test1 = X[train_index], X[test_index]
Y_train1, Y_test1 = Y[train_index], Y[test_index]
results = -1 * model_selection.cross_val_score(lm, X_train1,  Y_train1,scoring='neg_mean_squared_error', cv=kf)
print(results)
print(results.mean())
mse_test = mean_squared_error(X_test, Y_test)
print(mse_test)

#RMSE values interpretation test and train are similar no overfitting or underfitting the model performed well
print(math.sqrt(mse_test))
print(math.sqrt(results.mean()))


# Graph of the Training model
plt.scatter(X_train, Y_train, color = 'red')#plots scatter graph of COP COR  against PAUS for values in X_train and y_train
plt.plot(X_train, lm.predict(X_train), color = 'blue')#plots the graph of predicted PAUS against COP COR.
plt.title('SP000905974')
plt.xlabel('COP COR Quantity')
plt.ylabel('PAUS Quantity')
plt.show()#Show the graph

# Statistical Analysis of the training set with Statsmodels
X2 = sm.add_constant(X_train) # add a constant to the model
est = sm.OLS(Y_train, X2).fit()
print(est.summary()) # print the results

Answer 1

Assumptions

implement a simple linear regression to check the impact of variable A on B. I am not interested in predicting [...]

In one sense, linear regression produces predictions whether you want them or not. This has shades of meaning - possibly you don't care about interpolation/extrapolation/reconstruction, only metrics of correlation. Fine.

Is there anything fundamentally wrong with this approach to a simple Linear Regression?

This is not Cross Validated, the statistics sister site, so I won't go into depth about statistical validity, but I will show some simple metrics and a little bit of context.

I can't find your PL32_PMM_03_09_2018_SP_Level dataset anywhere public. Is PAUS the Physics of the Accelerating Universe Survey? Probably. Is COPCOR a reference to the CopCor coroplastics mold database referenced in Les Carnets de l'ACoSt 2016 (probably not), or the Cosmic ORigins Program and its COPAG (Cosmic Origins Program Analysis Group) - maybe? Answers to these questions are probably lost forever, so I have no choice but to demonstrate with a substitute dataset.

To make matters worse, your code would never have executed as posted, because (a) it refers to hard-coded filesystem paths - which shouldn't have been in there in the first place and will not be true for anyone else's computer - and (b) it is syntactically incorrect. You've removed all indentation.

I'm also not compelled to install sklearn or statsmodels, nor are they necessary for this analysis, so here is the set of assumptions I'll proceed with:

• Use an excerpt of the Open Asteroid Dataset.
• Skip all deserialisation (Excel, filesystem paths, etc.) concerns and load data inline
• Monovariate data, as you had.
• Only use methods from Numpy and Scipy, which are popular and well-documented
• Calling into linear regression methods is fine, and any "predictive" features will be ignored in favour of "descriptive" features

Demonstration

This will not be a course on linear algebra or descriptive statistics, only a demonstration that those statistics are easily accessible with what I consider more standard libraries than sklearn. Also: there are many, many more correlation tests you could run, but it's unlikely that they would be much more meaningful in context.

import numpy as np
import scipy.stats


def load_data() -> np.ndarray:
    """
    https://www.kaggle.com/datasets/basu369victor/prediction-of-asteroid-diameter
    columns H, rot_per
    excerpt of first 50 lines from CSV
    weak positive correlation
    """
    return np.array((
        (3.34,9.07417),
        (4.13,7.8132),
        (5.33,7.21),
        (3.2,5.34212766),
        (6.85,16.806),
        (5.71,7.2745),
        (5.51,7.139),
        (6.49,12.865),
        (6.28,5.079),
        (5.43,27.63),
        (6.55,13.7204),
        (7.24,8.6599),
        (6.74,7.045),
        (6.3,15.028),
        (5.28,6.083),
        (5.9,4.196),
        (7.76,12.27048),
        (6.51,11.57),
        (7.13,7.4432),
        (6.5,8.098),
        (7.35,8.1655),
        (6.45,4.1483),
        (6.95,12.312),
        (7.08,8.374),
        (7.83,9.9341),
        (7.4,13.11),
        (7,10.4082),
        (7.09,15.706),
        (5.85,5.3921),
        (7.57,13.686),
        (6.74,5.53),
        (7.56,9.448),
        (8.55,18.608),
        (8.51,12.15),
        (8.5,31.9),
        (8.46,9.93),
        (7.29,7.3335),
        (8.32,12.838),
        (6,5.138),
        (7,8.91),
        (7.12,5.988),
        (7.53,13.59),
        (7.93,5.76218),
        (7.03,6.422),
        (7.46,5.699),
        (8.36,21.04),
        (7.84,13.178),
        (6.9,11.89),
        (7.8,20.705),
    ))


def demo() -> None:
    data = load_data()
    n = len(data)
    x, y = data.T

    print(f'{n} samples loaded from Asteroid dataset')
    print()

    # a and b are the linear model fit, which we ignore.
    (a, b), (residual,), rank, (singular_1, singular_2) = np.linalg.lstsq(
        a=np.stack((x, np.ones_like(x)), axis=1),
        b=y,
    )
    print(f'''https://en.wikipedia.org/wiki/Rank_(linear_algebra)
Dimension of the vector space generated by the columns of the homogeneous linear fit matrix.
This must be 2 for a well-formed system.
rank = {rank}

https://en.wikipedia.org/wiki/Residual_sum_of_squares
Deviations predicted from actual empirical values of data. Measure of the discrepancy between the data and the estimation model.
SSR = {residual:.2f}

https://en.wikipedia.org/wiki/Singular_value
Positive, decreasing-order square root of the eigenvalues of the homogeneous linear fit matrix.
For a good intepretation, read https://sthalles.github.io/svd-for-regression/ .
singular values = {singular_1:.2f}, {singular_2:.2f}
''')

    # Identical to a, b above: ignore slope and intercept
    result = scipy.stats.linregress(x=x, y=y)

    print(f'''The following correlation tests come from
https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html

https://en.wikipedia.org/wiki/Pearson_correlation_coefficient
Measures linear correlation between two sets of data. Normalized measurement of the covariance, such that the result always has a value between −1 and 1.
Pearson r = {result.rvalue:.3f}

https://en.wikipedia.org/wiki/Wald_test
Assesses constraints on statistical parameters based on the weighted distance between the unrestricted estimate and its hypothesized value under the null hypothesis.
https://en.wikipedia.org/wiki/P-value
Probability of obtaining test results at least as extreme as the result actually observed, under the assumption that the null hypothesis is correct.
Wald p-value = {result.pvalue:.4f}

https://en.wikipedia.org/wiki/Standard_error
Standard deviation of the gradient's sampling distribution.
gradient stderr = {result.stderr:.3f}
Standard deviation of the intercept's sampling distribution.
intercept stderr = {result.intercept_stderr:.2f}''')


if __name__ == '__main__':
    demo()

49 samples loaded from Asteroid dataset

https://en.wikipedia.org/wiki/Rank_(linear_algebra)
Dimension of the vector space generated by the columns of the homogeneous linear fit matrix.
This must be 2 for a well-formed system.
rank = 2

https://en.wikipedia.org/wiki/Residual_sum_of_squares
Deviations predicted from actual empirical values of data. Measure of the discrepancy between the data and the estimation model.
SSR = 1387.97

https://en.wikipedia.org/wiki/Singular_value
Positive, decreasing-order square root of the eigenvalues of the homogeneous linear fit matrix.
For a good intepretation, read https://sthalles.github.io/svd-for-regression/ .
singular values = 48.88, 1.20

The following correlation tests come from
https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html

https://en.wikipedia.org/wiki/Pearson_correlation_coefficient
Measures linear correlation between two sets of data. Normalized measurement of the covariance, such that the result always has a value between −1 and 1.
Pearson r = 0.362

https://en.wikipedia.org/wiki/Wald_test
Assesses constraints on statistical parameters based on the weighted distance between the unrestricted estimate and its hypothesized value under the null hypothesis.
https://en.wikipedia.org/wiki/P-value
Probability of obtaining test results at least as extreme as the result actually observed, under the assumption that the null hypothesis is correct.
Wald p-value = 0.0105

https://en.wikipedia.org/wiki/Standard_error
Standard deviation of the gradient's sampling distribution.
gradient stderr = 0.649
Standard deviation of the intercept's sampling distribution.
intercept stderr = 4.49