I am playing with a toy data where the Simpson's paradox exists for two variables NO2 and temperature:

A scatter plot clearly shows that the correlation between NO2 and temperature was reversed when going from population level to subject level:

I am thinking of using lme to detect this, by specifying temperature as a fixed effect and for each subject a random intercept:

summary(lme(NO2~temperature, data=data.frame(data2), random = ~1| subject))


My expectation was that despite the negative correlation at the individual level, the fixed effect of temperature should reflect what is at the population level, i.e. a positive coefficient. However, the results showed the opposite as a negative coefficient for temperature:

My puzzle is:

1. It seems that the fixed effect estimation result captures the negative correlation between NO2 and temperature within each individual but why?

2. Is there a method to return both the correlations at the population level as well as the individual level? I tried to add a random slope for each subject like this:

summary(lme(NO2~temperature, data=data.frame(data2), random = ~0+temperature| subject)) And then it returned a positive coefficient. Why is that?

However, you can compare that model with another model without the individual factor, that is, NO2 ~ temperature, where coefficient for temperature will be positive. In fact, that is Simson's reversal: the sign of the coefficient of one predictor gets reversed when we take in account another (usually categorical) predictor.