4
$\begingroup$

I have measurements for daily space heating vs daily mean outdoor temperature for two different control strategies. The data is shown here:

enter image description here

I have also performed linear regression, of the form:

lm(Energy ~ Control * MeanOutdoorTemp)

This yields four coefficients:

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                170.8293     4.2083  40.594  < 2e-16 ***
ControlMPC                 -30.7044     4.9025  -6.263 1.38e-07 ***
MeanOutdoorTemp             -6.2924     1.6466  -3.821 0.000413 ***
ControlMPC:MeanOutdoorTemp   0.8211     1.7162   0.478 0.634709    

I'd like to test whether these two regression lines cross at zero energy. What would be the statistically correct way to do that?

$\endgroup$
3
  • $\begingroup$ A closely related question is discussed at stats.stackexchange.com/questions/15511/…. In the present case there are simpler solutions because you have specified beforehand the ordinate at which the crossing is to be tested. BTW, there's a good reason the lines shouldn't quite cross: daily fluctuations in temperature cause the true fit to be slightly curved for low heating energies. $\endgroup$
    – whuber
    Commented Jun 24, 2014 at 15:27
  • $\begingroup$ Physics says that wind speed and "characteristic length" of the house are also contributing factors. The Nusselt number drives the convective cooling coefficient, but is driven by the characteristic length (think radius of spherical house) nature of the fluid, flow speed etcetera. You are more likely to find a cubic (or other polynomial) relationship in house radius, and linear in flow speed, and linear in temperature difference. (en.wikipedia.org/wiki/…) $\endgroup$ Commented Jun 24, 2014 at 16:27
  • $\begingroup$ I notice you are trying to do both regressions at once using a dummy variable. Although that can work, in this case the scatterplot suggests the variances of the residuals depend on the control strategy. Consider, then, performing the regressions for each strategy separately. You still get four estimated coefficients (and they should have the same values as before) but now you will get two estimated residual variances to use in the test. Moreover, the estimates for one strategy are now independent of those for the other. $\endgroup$
    – whuber
    Commented Jun 24, 2014 at 16:40

1 Answer 1

2
$\begingroup$

The equation for the 1st line is $\beta_0 +\beta_1*x$ and part of the question is where this value equals 0, so you could set this equal to 0 and solve for $x$ to get $x=\frac{-\beta_0}{\beta_1}$. The equation for the 2nd line would then be $(\beta_0 + \beta_2) + (\beta_1 + \beta_3) * x$ and you can again set this to 0 and solve for $x$ to get $x=\frac{-(\beta_0 + \beta_2)}{(\beta_1 + \beta_3)}$. You can now test if the difference between these 2 is 0, or fit a confidence interval to the difference.

Unfortunately, working with ratios in cases like this can be very difficult (the distribution of a ratio of 2 normals can be a Cauchy if the denominator can be close to 0). One approach to get a reasonable approximation is to use simulation. This will assume that the estimates are normally distributed (standard assumption, but make sure that your residuals are normal enough or sample size large enough for this to be reasonable) You will need the covariance of your slope estimates (vcov(fit) in R where fit is the fitted model object). Now you can simulate slopes from a multivariate normal with mean equal to the estimated slopes and the given covariance matrix (in R use the mvrnorm function from the MASS package). The for each simulated set of values (4 slopes) compute the difference of the 2 equations above. Do this for a few thousand simulated values and then plot/summarize the results and see how they compare to 0.

Here is some sample R code to do a simulation:

x <- runif(50, 5, 10)
g <- rep(0:1, each=25)
y <- 0.2*x + 0.1*x*g + rnorm(50, 0, 0.1)
# y <- 0.2*x - 0.5*g + 0.1*x*g + rnorm(50, 0, 0.1)
plot(x,y, pch= g+1, xlim=c(0,10), ylim=range(0,y))

fit <- lm(y~x*g)
summary(fit)

tmp <- coef(fit)
abline( tmp[1], tmp[2], col='blue' )
abline( tmp[1] + tmp[3], tmp[2] + tmp[4], col='green' )

library(MASS)
sims <- mvrnorm(100000, tmp, vcov(fit))
diffs <- sims[,1]/sims[,2] - (sims[,1]+sims[,3])/(sims[,2]+sims[,4])
hist(diffs)
abline(v=0, col='red')
mean(diffs>0)
mean(diffs<0)
quantile(diffs, c(0.025, 0.975))

You can play around with this (the commented line gives different x intercepts by changing the y intercept of the 2nd group, run this line instead to see a difference, the uncommented line they do have the same x intercept (at 0)).

Try playing around with the parameters to see the effects. This can also give you a feel for power and what size of difference you can detect.

$\endgroup$
2
  • 1
    $\begingroup$ This looks like a good approach. I think you can do better by means of a little algebra to avoid the ratio problem. The test whether $y=\alpha_0+\alpha_1x$ and $y=\beta_0+\beta_1x$ intersect at $y=0$ comes down to whether $\alpha_0\beta_1-\alpha_1\beta_0=0$. Since $(\hat{\alpha_0},\hat{\alpha_1})$ and $(\hat{\beta_0},\hat{\beta_1})$ have independent bivariate normal distributions, we should be able to estimate the sampling distribution of $\hat{\alpha_0}\hat{\beta_1}-\hat{\alpha_1}\hat{\beta_0}$ without resorting to simulation. $\endgroup$
    – whuber
    Commented Jun 24, 2014 at 17:52
  • $\begingroup$ @whuber, good point. For testing equal to 0 we can get rid of the denominators. For estimating the difference with a confidence interval I am not sure if this will simplify things (but that is an extension of the original question). Products of normals are still more complicated than something that can be done with the general linear hypothesis. Someone with stronger theory skills than me should be able to work out a formal test, I still like the simulations (personal bias). $\endgroup$
    – Greg Snow
    Commented Jun 24, 2014 at 18:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.