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

Question

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)

Answer 1

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.

Answer 2

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)