<jupyter_start><jupyter_text>Where does model p-values come from?I've always wondered where p-values come from. I knew they would tell me how likely a model coefficient is telling me something about my data, but I never really knew how you get from point a to point b to point p-value. So here I've tried to describe the process for a linear model from nothing to coefficient p-value. I suppose this is basically true for any sort of GLM (like logistic regression, multi-level models, etc.). But in this case I'm doing it for a model with a single independent variable. Ordinary Least SquaresI'm not going to repeat all of wikipedia and every machine learning MOOC out there, and try to stick to the basics here. Ordinary least squares is a way to solve a linear model to find the best fitting model. Many modeling systems do this automatically (for e.g., [statsmodels](https://www.statsmodels.org/stable/index.html)). Basically, we want to take this equation:$$ \hat{y} = \hat{\beta} \hat{X} + b $$Where $\hat{y}$ is the dependent variable (the thing you want to predict) and $\hat{X}$ are all the variables (i.e., your data) you think might give you some prediction and perhaps explanation of where $\hat{y}$ comes from. Now this is just a model, it fits a plane to your data and suggests that plane explains something about the data or at least gives a prediction for what a given $y_i$ is given a bunch of $X_i$. In order to know what the shape of the plane is supposed to be, we have to know what the $\beta$s are. The way we figure out the best $\beta$s is by minimizing a loss function (i.e., least squares).$$ L(X) = \sum_{i}^{N} (y_i - p_i)^2 $$That is, given our actual y value in our data, how far away is our predicted model $\hat{y}$ and then do this for all the data being used to fit the model. Since this is the square error (basically) we want the error to be low because if its high we can't really trust what our model is telling us. So we have to minimize it. One way to minimize something is to take the derivative and set it equal to zero. This is typically done by [solving the normal equations](http://mathworld.wolfram.com/LeastSquaresFitting.html).So to bring this all down to earth, lets say we have some hypothesis that strong men's maximum deadlift is proportional to their total muscle mass. Let's also assume for argument's sake that we have some way to directly measure the muscle mass of a person and that our data set is controlled for things like height which may impact the total muscle mass.<jupyter_code>num_samples = 10000

weight_kg = np.random.normal(loc=150, scale=15, size=num_samples)
coef = 1.5
intercept = 200
max_dead_lift = intercept + coef * weight_kg + np.random.normal(loc=0, scale=15, size=num_samples)


fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))
ax1.scatter(weight_kg, max_dead_lift, alpha=0.5)
ax1.set_ylabel('maximum deadlift (kg)', fontsize=15)
ax1.set_xlabel('weight (kg)', fontsize=15)

ax2.hist(max_dead_lift, bins=100)
ax2.set_xlabel('max deadlift', fontsize=15)

ax3.hist(weight_kg, bins=100)
ax3.set_xlabel('weight (kg)', fontsize=15)
# note: ther aren't a lot of people who are >200kg, however eddie hall, the current record holder of deadlift (at 500kg) did weigh approximately 200kg when he won the record of a 500kg unassisted deadlift.<jupyter_output><empty_output><jupyter_text>So now that we have our fake data concerning strong men and their weights, we can try to fit a model between them iteratively. Given some starting value of the coefficient `beta` of our model:    max_deadlift ~ beta * weight + intercept    So basically we guess what the starting values of `beta` and `intercept` are, then we vary them by some step hoping that it gets better until it stops getting better. One way to guess is to set `beta` and `intercept` to zero and start increasing it (we can fudge a little knowing the relationship is positive from the graph). In this case, we probably have the computational power to explore the entire error space but if we had tons of data and this took a lot of time we might not be able to do do a for loop like I've done below and just run a while loop until we've minimized the mean square error enough.<jupyter_code>data_df = pd.DataFrame({'weight (kg)':weight_kg, 'max deadlift':max_dead_lift})
train_x, test_x, train_y, test_y = train_test_split(data_df['weight (kg)'], data_df['max deadlift'])
lr = LinearRegression()
lr.fit(train_x.to_frame(), train_y)
lr.score(test_x.to_frame(), test_y)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 5))

predicted = lr.predict(test_x.to_frame())
ax1.scatter(test_y, predicted, alpha=0.5, color='forestgreen')

residuals = test_y - predicted
ax2.scatter(test_y, residuals, alpha=0.5, color='forestgreen')

h = ax3.hist(residuals, bins=100, color='forestgreen')

ax1.set_xlabel('actual', fontsize=15)
ax1.set_ylabel('predicted', fontsize=15)
ax2.set_xlabel('actual', fontsize=15)
ax2.set_ylabel('residuals', fontsize=15)
ax3.set_xlabel('residuals', fontsize=15)<jupyter_output><empty_output><jupyter_text>So basically, we have a linear model now that tells us this relationship of weight to max deadlift is probably telling us something about something. However, in order to have more believability of our model, we can perform various tests on the coeficients themselves. For example, we can perform a two tailed t-test between a model where the coefficient in question is 0 or the fitted value **[[citation](http://www.eco.uc3m.es/~ricmora/miccua/materials/DavidsonMacKinnon_2004_Ch4.pdf)]**. We will do this below.<jupyter_code>metrics.r2_score(test_y.values.reshape(-1, 1),lr.predict(test_x.to_frame()))
lr.coef_[0]
lr.intercept_
coefs = np.linspace(lr.coef_[0]-lr.coef_[0], lr.coef_[0]+lr.coef_[0], 100)
intercepts = np.linspace(lr.intercept_-lr.intercept_, lr.intercept_+lr.intercept_, 100)
y_pred = []
for m, b in itertools.product(coefs, intercepts):
    y = m * test_x + b
    y_pred.append((m, b, y))
    
r2s = []
for m, b, y in y_pred:
    r2 = metrics.r2_score(test_y.values.reshape(-1, 1), y)
#     ttest = stats.ttest_ind(y, yzero)
    ttest = stats.ttest_1samp(a=[b, m], popmean=0)
    r2s.append(np.array([m, b, r2, ttest.statistic, ttest.pvalue]))

# ttests = [stats.ttest_ind(y[2], yzero).pvalue for y in y_pred]
# ttests = [stats.ttest_1samp(a=[b,m], popmean=0).pvalue for y in y_pred]
r2_df = pd.DataFrame(np.array(r2s), columns=['coef', 'intercept', 'r2', 'ttest', 't-pval'])
r2_df['log_r2'] = r2_df.r2.apply(lambda r: np.log10(r))
r2_space = r2_df.pivot_table(index='coef', columns='intercept', values='r2')
r2_df.describe()
fig, ax = plt.subplots(figsize=(10, 9))
cbar = ax.pcolormesh(r2_space.columns, r2_space.index, r2_space.values, cmap='nipy_spectral', vmin=-0.7, vmax=0.7)
fig.colorbar(cbar)
ax.set_ylabel('Coefficient', fontsize=15)
ax.set_xlabel('Intercept', fontsize=15)
fig.get_axes()[1].set_ylabel('R2', fontsize=35)
p_space = r2_df.pivot_table(index='coef', columns='intercept', values='t-pval')

fig, ax = plt.subplots(figsize=(10, 9))
cbar = ax.pcolormesh(p_space.columns, p_space.index, p_space.values, cmap='gnuplot2_r')
fig.colorbar(cbar)
ax.set_ylabel('Coefficient', fontsize=15)
ax.set_xlabel('Intercept', fontsize=15)
fig.get_axes()[1].set_ylabel('p-value', fontsize=35)
yzero = 0 * test_x + lr.intercept_

# ttests = [t for t in stats.ttest_ind(y, yzero) for y in ypred]
for y in y_pred:
    stats.ttest_ind(y[2], yzero)
    
ttests = [stats.ttest_ind(y[2], yzero).pvalue for y in y_pred]
len(ttests)
ypred = lr.coef_[0]*test_x + lr.intercept_
stats.ttest_1samp(a=[lr.coef_[0], lr.intercept_], popmean=0)
lr.intercept_
import statsmodels.api as sm
sm_train_x = sm.add_constant(train_x)
model = sm.OLS(train_y, sm_train_x)
results = model.fit()
results.summary()<jupyter_output>/home/john/anaconda3/envs/research/lib/python3.6/site-packages/numpy/core/fromnumeric.py:2389: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
  return ptp(axis=axis, out=out, **kwargs)
{'fraction_non_alphanumeric': 0.11195611519760629, 'fraction_numerical': 0.023064455803515772, 'mean_word_length': 5.147126436781609, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 13, 'https://': 1, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '531468', 'license': 'MIT', 'n_tokens_mistral': 2828, 'n_tokens_neox': 2693, 'n_words': 973, 'path': 'ipython-notebooks/where do OLS p-values come from.ipynb', 'repo_name': 'mnky9800n/python-plotting', 'size': '8021'}