# bias of p-value analysis in a not well fitting model

In the report http://www.stat.berkeley.edu/~breiman/wald2002-3.pdf Breiman says:

Three decades ago many statisticians and quantitative social scientists were enamored of multilinear regression and its theory of hypothesis testing on the coefficients. Every statistical package had a regression program variable selection program based on F-to delete and F to enter. It was almost impossible to get a paper published unless you showed that a certain coefficient was signifigant at the 5% level. This was regardless of how well the linear model fit the data and little effort was made to find out. Many conclusions were undoubtedly wrong, and I don't think statisticians now-a-days dispute the error of these ways

I would like to have an example obtained with simulated data where a linear regression show that a coefficient is significant while it is not true. And I would like that the error will arise from the fact that the linear model is not fitting well the data.

EDIT: R or python code is welcome

EDIT: The example should make use of at least 1000 samples

EDIT: I found a possible example with simulated data.

      n = 1000
x1 = seq(-10, 10, length.out = 1000)
x2 = sin(x1)/x1 + rnorm(n = n, mean = 0, sd = 2)
y = sin(x1)/x1 + rnorm(n = n, mean = 0, sd = .1)
lr = glm(y~x1+x2) #wrong model
summary(lr)
z = sin(x1)/x1
lr = glm(y~z+x2) #real model
summary(lr)

plot(x1, y, col='blue')
points(x1, x2)
points(x1, z)

• "The example should make use of at least 1000 samples"... so is this homework? Otherwise, why 1000 as a specific lower limit? Mar 4, 2015 at 15:12
• @Glen_b not it is not homework. It is a way to avoid example that are not if interest for me Mar 4, 2015 at 15:13
• Why would 1000 be of interest, but 900 not be of interest? Mar 4, 2015 at 15:14
• The paragraph starting with I would like to... is not completely clear, judging from the discussion in the comments below my answer. Maybe you could specify more exactly what you are looking for? Mar 4, 2015 at 15:14
• The logic of this request is mysterious, because it is problematic to ascribe a truth value to the "significance" of a coefficient in a linear regression which is understood to be an incorrect model. Suppose the linear regression indicates a coefficient is significant. That could--and often does--happen because the corresponding variable is capturing, albeit imperfectly, some real form of variation in the data. We should not be so hasty to equate "imperfect" and "not true," however.
– whuber
Mar 4, 2015 at 17:51

Try a regression of $y=e^{2x}$ on $x$ for $x$ distributed uniformly in the interval $[0,1]$.

In R it can be coded as follows:

set.seed(1)
x=runif(10^3)
y=exp(2*x)
lm1=lm(y~x)
summary(lm1)


The true model is nonlinear, the estimated model is linear but the coefficient is highly significant.

Update: An example of spurious correlation between integrated processes:

set.seed(1); y=cumsum(rnorm(10^3))
set.seed(2); x=cumsum(rnorm(10^3))
lm1=lm(y~x)
summary(lm1)


The two variables $y$ and $x$ are unrelated but there is a spurious correlation between them, and the regression coefficient is highly significant. Many real world examples can be found here.

• Thanks for answering. In this case the linear model is able to correctly classify x as a relevant variable even if the relation of x and y is not linear. I would like a situation where x seems relevant while it is not or when it does not seem relevant while it is. Mar 4, 2015 at 14:40
• Then spurious correlations from unrelated integrated time series is the answer. Here are bunches of examples where completely irrelevant regressors would be statistically significant. Mar 4, 2015 at 14:44
• I am familiar with these examples. The spurious correlation there is caused mainly by the lack of data. While I am looking for error due to the use of the wrong model. Mar 4, 2015 at 14:46
• I am afraid spurious correlations between integrated processes are not due to lack of data. t-statistics have different distributions there so the phenomenon remains even in the asymptotics. Mar 4, 2015 at 14:54
• I am curious about why someone down voted this answer.. Mar 4, 2015 at 15:18

I found a possible example but it is a bit artificial and maybe it would be difficult to find something similar in real life

suppose $y = e^{x_1}+err$ and $x_2 =e^{x_1}+ERR$ where $ERR$ and $err$ are two measurement errors and $ERR>>err$. In this case the model $y=x1+x_2$ will consider $x_2$ as relevant variable while the model $y=e^{x_1}+x_2$ will consider $x_1$ as relevant.

> n = 1000
> x1 = seq(-10, 10, length.out = 1000)
> x2 = exp(x1) + rnorm(n = n, mean = 0, sd = 4)
> y = exp(x1) + rnorm(n = n, mean = 0, sd = 1) #y is a function of x1 only
> lr = glm(y~x1+x2) #the wrong model is used
> summary(lr)

Call:
glm(formula = y ~ x1 + x2)

Deviance Residuals:
Min        1Q    Median        3Q       Max
-11.9787   -2.8482    0.0355    2.7221   11.9990

Coefficients:
Estimate Std. Error   t value Pr(>|t|)
(Intercept) -2.347e-01  1.383e-01    -1.697    0.090 .
x1          -1.056e-02  2.610e-02    -0.404    0.686
x2           1.000e+00  4.521e-05 22116.688   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 16.60931)

Null deviance: 1.1126e+10  on 999  degrees of freedom
Residual deviance: 1.6559e+04  on 997  degrees of freedom
AIC: 5652.8

Number of Fisher Scoring iterations: 2

> lr = glm(y~exp(x1)+x2)
> summary(lr)

Call:
glm(formula = y ~ exp(x1) + x2) #the real model is used

Deviance Residuals:
Min       1Q   Median       3Q      Max
-3.2881  -0.6651   0.0234   0.6947   3.1845

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.003736   0.033282  -0.112    0.911
exp(x1)      0.990101   0.007912 125.147   <2e-16 ***
x2           0.009886   0.007911   1.250    0.212
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.99421)

Null deviance: 1.1126e+10  on 999  degrees of freedom
Residual deviance: 9.9123e+02  on 997  degrees of freedom
AIC: 2837.1

Number of Fisher Scoring iterations: 2

>


EDIT:

A very similar but possibly better example

n = 1000
x1 = seq(-10, 10, length.out = 1000)
x2 = sin(x1)/x1 + rnorm(n = n, mean = 0, sd = 2)
y = sin(x1)/x1 + rnorm(n = n, mean = 0, sd = .1)
lr = glm(y~x1+x2) #wrong model
summary(lr)
z = sin(x1)/x1
lr = glm(y~z+x2) #real model
summary(lr)

plot(x1, y, col='blue')
points(x1, x2)
points(x1, z)

• This example does not satisfy the requirement in the original post. Finding a significant $x_2$ in the regression $y=\beta_0+\beta_1 x_1+\beta_2 x_2+\epsilon$ makes perfect sense as $y=0+0 \cdot x_1+1 \cdot x_2+\epsilon$ where $\epsilon=err-ERR$. Thus you actually have a correct model with one extra irrelevant regressor which is found to be insignificant just as it should be. This violates coefficient is significant while it is not true if I decipher the it is not true correctly. However, it is not true should be treated with care (see e.g. a valid comment by whuber under the OP). Mar 4, 2015 at 21:19