Physically implausible results in linear regression with colinearities

Question

While developing a model (a Poisson regression, but this is not the topic of this post), I stumbled upon a physically implausible relationship between some variables.

I have ground temperature data for a European country. Nothing exotic:

I have also collected the Moon brightness (APmag) and Sun-Earth-Moon angle (SOT.angle) data from NASA. These follow an obvious pattern and are closely, but not exactly linearly related:

When I run a linear regression on these data, I find no relationship between the temperature and the brightness, angle, or both ($p > 0.8$, $R^2 \approx 0$). Again, no surprise here. However, if I take the interaction into account I suddenly obtain highly "significant" relationship ($p < 10^{-11}$, $R^2 \approx 0.14$):

Call:
lm(formula = TG ~ SOT.angle * APmag, data = hd2018)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.3948  -5.5808   0.5258   5.8828  16.9561 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     72.54770    8.73268   8.308 2.04e-15 ***
SOT.angle        2.46045    0.31466   7.819 6.00e-14 ***
APmag           15.72812    2.11922   7.422 8.50e-13 ***
SOT.angle:APmag  0.13507    0.01722   7.846 5.01e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.47 on 358 degrees of freedom
Multiple R-squared:  0.1468,    Adjusted R-squared:  0.1396 
F-statistic: 20.53 on 3 and 358 DF,  p-value: 2.692e-12

And that was just with the data for one year. If I take 20 years, I get p-values in the range $10^{-100}$.

In medicine, where I do my research, this would be considered proof beyond any doubt, but physically, it is obvious that the Moon cannot influence ground temperature on Earth*. My guess is it has to do something with the near-linear relationship between the brightness and the angle, but I can't say I really understand the mechanism.

In this case, I was lucky that my domain knowledge sufficed to identify this as a false relationship, but I fear there might be other false relationships for which we lack such domain knowledge.

Is there a statistical approach to avoid falling in such traps? I'd appreciate any clarification: mathematical, graphical, intuitive...

* Well, yes, maybe in the range $\ll 10^{-6} K$, which is far below the measurement precision and the noise level. I'm pretty sure the above analysis didn't discover that effect.

Update:

The actual temperature has nothing to do with the relationship. I get "significant" results ($p$ ranging from $10^{-3}$ to $10^{-13}$) for all of the following artificial "temperature" curves:

Answer 1

There is a confounder: The Earth's position on its orbit around the Sun determines on one hand the season (and, consequently, the ground temperature) and, on the other hand, the amount of light the Moon is receiving. There are subtle traces in the data.

The Sun-Earth-Moon angle and the Moon brightness are tightly related, which is not surprising. The Moon brightness, as observed from the Earth, is a direct consequence of the illuminated area we can observe. That area is given by the above angle. We can draw a smooth curve (3rd degree polynomial) through the points:

There is some "noise" (residuals) around the curve, but its level is so low that we might be tempted to ignore it. That, however, would be a big mistake! We are talking celestial mechanics here, not social sciences. There should be no noise at all! Everything is deterministic, and must have an explanation.

If we plot the noise against the date, another pattern emerges:

In summer (in the northern hemisphere), the residuals tend to be larger (i.e. the true Moon brightness is lower than the modeled one). Most probably, this is a result of the Earth's orbit around the Sun being elliptical. During summer, the Earth+Moon-system is farther away from the Sun, and the Moon gets somewhat less light (the Earth too, but this is negligible here).

There is also an oscillatory pattern (13 periods a year), probably due to the Moon's orbit around the Earth (also elliptic, and inclined towards the Earth's orbit), but that one doesn't account for the observed effect. This can be shown by creating an artificial "Moon brightness", modeled as a quadratic function of the date (the blue parabola in the above figure) and a cubic polynomial of SOT.angle (the smooth line in the first figure). This artificial "Moon brightness" still has a "significant" effect on the ground temperature TG!

In fact, the ground temperature can be predicted from the Moon brightness residual alone:

Call:
lm(formula = TG ~ AP.resid, data = hd2)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.9810  -5.5160   0.3846   6.4040  16.7107 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  11.6559     0.3987  29.233  < 2e-16 ***
AP.resid     29.0386     4.4487   6.527 2.26e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.618 on 363 degrees of freedom
Multiple R-squared:  0.105, Adjusted R-squared:  0.1026 
F-statistic: 42.61 on 1 and 363 DF,  p-value: 2.256e-10

It is not a time-series issue. In this case, of course, the positions of the celestial bodies are functions of time, but this is not the cause for the relationship in question. We could construct a laboratory experiment, with a light source, spherical mirrors, and light detectors, and perform measurements in any order we desire, and we'd still obtain the same relationship.

Answer 2

I suspect that a major problem is the time-series nature of the data.

As the Wikipedia page says, inference with linear regression assumes that "errors of the response variables [around the values predicted by the linear model] are uncorrelated with each other." That is often not the case with time series, for which the errors around the predicted values for observations nearby in time are often correlated with each other. This autocorrelation means that the independence assumptions underlying interpretation of p-values no longer hold. While that might not have affected your model without the interaction term, including an interaction term between 2 predictors that are already highly correlated might have led to this problem.

Another thing to think about in models with interaction terms (even absent autocorrelation issues) is that the p-values for coefficients other than the interaction term often aren't very useful. As usually reported, they represent coefficients for a situation in which all categorical predictors are at their reference levels and continuous-valued predictors have values of 0. For example, the coefficient and p-value reported for SOT.angle is for a situation in which APmag is 0, which is far out of the range of actual APmag values. So do you really care whether the coefficient for SOT.angle is different from 0 under those circumstances? That's all that the p-value for that coefficient tells you. That doesn't explain the low p-values for the model as a whole or for the interaction term itself; those probably come from autocorrelations among error terms in your model.