import os
import statsmodels.api as sm
if not os.getenv("GEMINI_API_KEY") and not os.getenv("GOOGLE_API_KEY"):
raise EnvironmentError("Please set your GEMINI_API_KEY or GOOGLE_API_KEY environment variable.")
# 1. Load data and fit a model (Duncan's occupational prestige data)
duncan_data = sm.datasets.get_rdataset("Duncan", "carData")
y = duncan_data.data["prestige"]
X = sm.add_constant(duncan_data.data[["income", "education"]])
model = sm.OLS(y, X).fit()Get started (Python)
Installation
From PyPI (once published):
pip install statlingoOr install the development version directly from GitHub:
pip install "git+https://github.com/bgreenwell/statlingo.git#subdirectory=python"You’ll also need a chatlas-supported provider SDK (e.g. openai for ChatOpenAI) and statsmodels for fitting the example model below.
Quick example
First, import the required packages, check environment variables, load data, and fit an OLS model:
This checks for local API credentials and fits an ordinary least squares regression model of occupational prestige on income and education using statsmodels.
Next, initialize the chatlas chat client using Google Gemini:
from chatlas import ChatGoogle
# 2. Create a chatlas Chat client using Google Gemini (never mutated by explain())
client = ChatGoogle()This initializes a Google Gemini chat client from chatlas. This client is deep-copied internally within explain() to prevent mutating your active chat session.
Generate the natural language explanation using explain(), targeting a student audience:
from statlingo import explain
# 3. Get a high-level explanation of the model results
explanation = explain(model_object=model, client=client, audience="student")
print(explanation["text"])This output represents the results of an Ordinary Least Squares (OLS) Linear Regression model. This type of model aims to find a linear relationship between a continuous response variable and one or more predictor variables. The goal is to predict the average value of the response variable based on the values of the predictors.
## Model Appropriateness
Based solely on this output, a linear regression model *could* be appropriate if your response variable (`prestige`) is continuous and you hypothesize a linear relationship with your predictors (`income` and `education`). However, without more context about your data, specific research question, and how the variables are measured, it's impossible to definitively confirm the model's appropriateness or fully assess if all underlying assumptions are met.
## Interpretation of the Output
Let's break down the different sections of this output.
### Model Specification and Overall Fit
This top section provides information about how the model was set up and some overall measures of its performance.
* **Dep. Variable: prestige**: This tells us that `prestige` is your *dependent variable* (also called the response variable). This is the outcome you are trying to predict or explain.
* **Model: OLS**: Confirms that an Ordinary Least Squares method was used, which is standard for linear regression.
* **Method: Least Squares**: This is the algorithm used to find the "best-fitting" line by minimizing the sum of the squared differences between the observed and predicted values.
* **No. Observations: 45**: You have 45 data points (or individuals/cases) in your dataset used for this analysis.
* **Df Residuals: 42**: The *degrees of freedom for the residuals*. This is calculated as `No. Observations - No. of Coefficients Estimated`. In this case, $45 - 3 = 42$ (3 coefficients: `const`, `income`, `education`).
* **Df Model: 2**: The *degrees of freedom for the model*. This is the number of predictor variables in your model (excluding the intercept), which are `income` and `education`.
* **R-squared: 0.828**: This is the *coefficient of determination*. It tells you that approximately 82.8% of the variance (or variability) in `prestige` can be explained by your predictor variables (`income` and `education`) in this model. This is a relatively high R-squared, suggesting the model explains a large portion of the variation in prestige.
* **Adj. R-squared: 0.820**: The *Adjusted R-squared* is a modified version of R-squared that accounts for the number of predictors in the model. It tends to be a more reliable measure when comparing models with different numbers of predictors, as it penalizes models for including unnecessary predictors. It's slightly lower than R-squared, which is typical.
* **F-statistic: 101.2**: This is the F-statistic for the overall model. It tests the null hypothesis that *all* of your regression coefficients (excluding the intercept) are simultaneously equal to zero. In simpler terms, it tests if your model with predictors is significantly better at explaining `prestige` than a model with no predictors (just an intercept).
* **Prob (F-statistic): 8.65e-17**: This is the p-value associated with the F-statistic. A very small p-value (like `8.65e-17`, which is essentially 0) indicates strong evidence to reject the null hypothesis. This means that your model, as a whole, is statistically significant, and at least one of your predictor variables (`income` or `education`) is significantly related to `prestige`.
* **Log-Likelihood, AIC, BIC**: These are criteria used for model comparison, especially when comparing non-nested models or models fit with different methods. Lower values generally indicate a better fit, but they are most useful in a comparative context.
### Coefficients Table
This is the core of your model's results, showing the estimated effect of each predictor on the response variable.
| Term | coef | std err | t | P>|t| | [0.025 | 0.975] |
| :---------- | :------ | :------ | :------ | :---- | :----- | :----- |
| `const` | -6.0647 | 4.272 | -1.420 | 0.163 | -14.686| 2.556 |
| `income` | 0.5987 | 0.120 | 5.003 | 0.000 | 0.357 | 0.840 |
| `education` | 0.5458 | 0.098 | 5.555 | 0.000 | 0.348 | 0.744 |
* **`const` (Intercept)**:
* **coef: -6.0647**: This is the estimated *intercept* of the regression line. It represents the predicted value of `prestige` when `income` and `education` are both zero. In many real-world scenarios, a zero value for predictors like income or education might not be meaningful or even possible. Therefore, interpreting the intercept directly might not always be practical or relevant, but it's crucial for defining the regression line.
* **P>|t|: 0.163**: With a p-value of 0.163, which is greater than common significance levels (e.g., 0.05), we do not have sufficient evidence to conclude that the true intercept is statistically different from zero.
* **`income`**:
* **coef: 0.5987**: This is the estimated *regression coefficient* for `income`. It means that, *holding `education` constant*, for every one-unit increase in `income`, the predicted `prestige` increases by approximately 0.5987 units. (If you know the units of `income` and `prestige`, you can make this interpretation more specific).
* **std err: 0.120**: This is the *standard error* of the `income` coefficient. It measures the precision of the `income` estimate. A smaller standard error indicates a more precise estimate.
* **t: 5.003**: This is the *t-statistic* for `income`. It's calculated by dividing the coefficient estimate by its standard error ($0.5987 / 0.120 \approx 5.003$). This statistic tests the null hypothesis that the true coefficient for `income` is zero (meaning `income` has no linear relationship with `prestige`).
* **P>|t|: 0.000**: This is the p-value associated with the `income` t-statistic. A p-value of 0.000 (which is less than 0.05) indicates very strong evidence against the null hypothesis. Therefore, `income` is a statistically significant predictor of `prestige` in this model.
* **[0.025 0.975]: [0.357 0.840]**: This is the 95% confidence interval for the `income` coefficient. It means we are 95% confident that the true effect of a one-unit increase in `income` on `prestige` (holding `education` constant) lies between 0.357 and 0.840 units. Since this interval does not include zero, it reinforces that `income` is a statistically significant predictor.
* **`education`**:
* **coef: 0.5458**: This is the estimated regression coefficient for `education`. It means that, *holding `income` constant*, for every one-unit increase in `education`, the predicted `prestige` increases by approximately 0.5458 units.
* **std err: 0.098**: The standard error of the `education` coefficient.
* **t: 5.555**: The t-statistic for `education`.
* **P>|t|: 0.000**: The p-value for `education`. A p-value of 0.000 (less than 0.05) indicates strong evidence against the null hypothesis, suggesting that `education` is also a statistically significant predictor of `prestige`.
* **[0.025 0.975]: [0.348 0.744]**: This is the 95% confidence interval for the `education` coefficient. We are 95% confident that the true effect of a one-unit increase in `education` on `prestige` (holding `income` constant) lies between 0.348 and 0.744 units. This interval also does not include zero, confirming its statistical significance.
### Bottom Diagnostics Block
This section provides additional diagnostic tests to help assess whether the model's assumptions are being met.
* **Omnibus: 1.279, Prob(Omnibus): 0.528**: The Omnibus test (often D'Agostino-Pearson's K-squared test) tests the normality of the residuals. A high p-value (0.528) suggests that we *do not* reject the null hypothesis of normally distributed residuals. This is a good sign for the normality assumption.
* **Skew: 0.155**: This measures the asymmetry of the residual distribution. A value close to zero indicates symmetry. 0.155 is close to zero, suggesting little skewness.
* **Kurtosis: 3.426**: This measures the "tailedness" of the residual distribution (how many outliers there are). For a normal distribution, kurtosis is 3. A value of 3.426 is reasonably close to 3, suggesting the tails are not excessively heavy or light compared to a normal distribution.
* **Durbin-Watson: 1.458**: This statistic tests for *autocorrelation* in the residuals (i.e., whether the residuals are independent of each other). A value close to 2 indicates no autocorrelation. Values below 1 or above 3 typically suggest issues. A value of 1.458 suggests there might be some positive autocorrelation, but it's not extremely strong. This assumption is particularly important for time series data.
* **Jarque-Bera (JB): 0.520, Prob(JB): 0.771**: Another test for the normality of residuals, combining skewness and kurtosis. A high p-value (0.771) suggests that we *do not* reject the null hypothesis of normally distributed residuals. This reinforces the finding from the Omnibus test.
* **Cond. No.: 163**: The *Condition Number* assesses multicollinearity, which is when predictor variables are highly correlated with each other. A condition number above 30 indicates moderate to severe multicollinearity. Your value of 163 suggests that there might be a significant multicollinearity issue between `income` and `education`. This can make the individual coefficient estimates less stable and harder to interpret, even if the overall model (F-statistic) is significant.
## Suggestions for Checking Assumptions
While some diagnostic statistics are provided, it's highly recommended to use graphical methods to visually inspect the model's assumptions. These plots often reveal more nuanced issues than statistical tests alone.
1. **Residuals vs. Fitted Values Plot**:
* **Purpose**: To check for *linearity* and *homoscedasticity* (constant variance of errors).
* **What to look for**: A random scatter of points around the horizontal line at zero, with no discernible pattern (like a curve or a cone shape).
* **Issues**: A clear pattern (e.g., a U-shape) suggests non-linearity. A fanning-out or fanning-in pattern suggests heteroscedasticity (non-constant variance).
2. **Normal Q-Q Plot of Residuals**:
* **Purpose**: To check the *normality* of the residuals.
* **What to look for**: Points that closely follow the diagonal line.
* **Issues**: Points deviating significantly from the line, especially at the tails, suggest non-normality.
3. **Scale-Location Plot (Square Root of Standardized Residuals vs. Fitted Values)**:
* **Purpose**: Another way to check for *homoscedasticity*.
* **What to look for**: A horizontal line with randomly scattered points.
* **Issues**: A clear pattern (e.g., a rising or falling trend) indicates heteroscedasticity.
4. **Residuals vs. Leverage Plot**:
* **Purpose**: To identify *influential observations* that might unduly affect the regression results.
* **What to look for**: Points far from the main cluster, especially those outside Cook's distance contours (if plotted).
Given the `Cond. No.` of 163, you should also investigate multicollinearity further, perhaps by calculating the Variance Inflation Factors (VIFs) for your predictors. If VIF values are high (e.g., above 5 or 10), it confirms a multicollinearity problem. For the Durbin-Watson statistic, if your data has a time component or a natural ordering, you would particularly want to investigate the potential autocorrelation.
***
*Caution: This explanation was generated by a Large Language Model. Users should critically review the output and consult additional statistical resources or experts to ensure correctness and a full understanding.*
This generates a dictionary output whose text key contains the formatted student-oriented OLS explanation.
Finally, request diagnostic code recommendations based on the generated explanation using suggest_code():
from statlingo import suggest_code
# 4. Suggest next diagnostic code steps based on the explanation
print(suggest_code(explanation))# Next Steps: Suggested Python Coding Diagnostics
# 1. Plot Residuals vs Fitted values
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
# Assuming 'model' is your fitted OLSResults object
fig, ax = plt.subplots(figsize=(8, 5))
sns.residplot(x=model.fittedvalues, y=model.resid, lowess=True,
scatter_kws={'alpha': 0.5},
line_kws={'color': 'red', 'lw': 1, 'alpha': 0.8}, ax=ax)
ax.set_title('Residuals vs Fitted')
ax.set_xlabel('Fitted values')
ax.set_ylabel('Residuals')
plt.show()
# 2. Test for Homoscedasticity (Breusch-Pagan)
from statsmodels.stats.diagnostic import het_breuschpagan
bp_test = het_breuschpagan(model.resid, model.model.exog)
labels = ['Lagrange multiplier statistic', 'p-value', 'f-value', 'f p-value']
print(dict(zip(labels, bp_test)))
# 3. Test for Autocorrelation (Durbin-Watson)
from statsmodels.stats.stattools import durbin_watson
dw = durbin_watson(model.resid)
print(f'Durbin-Watson statistic: {dw}')
# 4. Check Multicollinearity (VIF)
from statsmodels.stats.outliers_influence import variance_inflation_factor
for i in range(1, model.model.exog.shape[1]):
vif = variance_inflation_factor(model.model.exog, i)
print(f'VIF for predictor {i}: {vif}')
***
*Caution: These code suggestions were generated by a Large Language Model. Users should critically review the code to ensure correctness and applicability before executing.*
This prints recommended Python code using matplotlib, seaborn, and statsmodels.stats to check homoscedasticity, autocorrelation, and multicollinearity.