# Confidence intervals for mean response do not agree using predict.lm() in R

I have $$n$$ data points $$y_i,x_i$$ and the model $$y_i=\alpha x_i^\beta + \epsilon_i$$ where $$\log{\epsilon_i} \in N(0,\sigma^2)$$. Thus, I define $$z_i=\log{y_i}$$ and $$w_i=\log{x_i}$$ and my normal linear model becomes $$z_i=\log{\alpha}+\beta w_i+\log{\epsilon_i}$$. I'd like to construct a 99% confidence interval for $$E(\hat{y}(x^0))$$, i.e. the mean response given some value $$x_0$$.

For a model of the form,$$y_i=\theta_0+\theta_1 (x_{1i}-\bar{x}_{1\cdot})+...+\theta_p (x_{pi}-\bar{x}_{p\cdot}) + \epsilon_i \ , \tag{1}$$ where $$\bar{x}_{j\cdot}=n^{-1}\sum_{i=1}^n x_{ij}$$, the confidence interval at some point $$x^0=(x_1^0,...,x_p^0)$$ is given by,

$$\hat{y}(x^0)\pm t_{\alpha/2}(n-p-1)s_p\sqrt{\frac{1}{n}+x_v^T(X^TX)^{-1}x_v} \ , \tag{2}$$

where $$s_p$$ is the biased corrected MLE of $$\hat{\sigma}$$, that is

$$s_p=\frac{Q(\hat{\theta})}{n-p-1} \ , \tag{3}$$

and $$x_v=(x_0^1-\bar{x}_{1\cdot},...,x_0^p-\bar{x}_{p\cdot})$$.

In Python I'd do the following (with $$x_0=6$$ and $$t_{0.005}(8)\approx 3.355$$):

import numpy as np

x=np.array([2,3,4,5,6,7,8,9,10,11])                 #n=10, p=1
y=np.array([2.1,4,3.7,4.5,5,4.8,5.1,5.7,5.7,5.6])

z=np.log(y)
w=np.log(x)

theta_0=np.mean(z)                                  #\theta_0 in equation (1) (=\log{\alpha}+\hat{\beta}*\bar{w})

num=z*(w-np.mean(w))
den=(w-np.mean(w))**2
theta_1=np.sum(num)/np.sum(den)                     #\theta_1 in equation (1) (=\beta)

Q=np.sum( (z-(theta_0+theta_1*(w-np.mean(w))))**2 )
s_p=np.sqrt(Q/8)

z_log6=(theta_0-theta_1*np.mean(w))+theta_1*np.log(6)
lhs=z_log6 - 3.355*s_p*np.sqrt( 1/10 + 1/np.sum((w-np.mean(w))**2) ) #equation (2)
rhs=z_log6 + 3.355*s_p*np.sqrt( 1/10 + 1/np.sum((w-np.mean(w))**2) ) #equation (2)
print(lhs,rhs)
print(np.exp(lhs),np.exp(rhs))


With output:

1.2313777318582835 1.8006635458783546
3.4259463217029915 6.053663015157194


In R and using predict.lm():

x<-c(2,3,4,5,6,7,8,9,10,11)
y<-c(2.1,4,3.7,4.5,5,4.8,5.1,5.7,5.7,5.6)

w<-log(x)
z<-log(y)

fit<-lm(z~w)
summary(fit)

newdata<-data.frame(w=log(6))
print(predict.lm(fit,newdata,interval = "confidence", level = 0.99))


With output:

       fit      lwr      upr
1 1.516021 1.383078 1.648963


The first confidence intervals for both outputs, i.e. for $$E(\hat{z}(\log{6}))$$ do not agree (and hence not the second either. Why is that?

• Looking at www2.stat.duke.edu/~tjl13/s101/slides/unit6lec3H.pdf slide 5, I think there should be a $(n-1)$ in the denominator. Also, check that you using the same statistical significance $a$. – Fiodor1234 Feb 1 at 19:17
• Thanks for the reply. I think it is correct. In the Python code, I have in the denominator np.sum((w-np.mean(w))**2) ). The denominator is the unbiased sample variance less than a factor $(n-1)$, so I could replace it with the unbiased sample variance ($s_x^2$ in the slides) but then add a factor $(n-1)$ as well. – schn Feb 1 at 19:33

## 1 Answer

There's an error in your formula in the Python program.
The formula in R is shown here:

> print(predict.lm(fit,newdata,interval = "conf", level = 0.99))
fit      lwr      upr
1 1.516021 1.383078 1.648963
> 1.516021+c(-1,1)*qt(0.995,8)*0.1249*sqrt((1/10+(log(6)-mean(w))^2/sum((w-mean(w))^2)))
[1] 1.383081 1.648961


In your Python code, you should replace the 1 in the numerator of the second term in the sqrt with the right thing from above.

lhs=z_log6 - 3.355*s_p*np.sqrt( 1/10 + (log(6)-np.mean(w))**2/np.sum((w-np.mean(w))**2) ) #equation (2)
rhs=z_log6 + 3.355*s_p*np.sqrt( 1/10 + (log(6)-np.mean(w))**2/np.sum((w-np.mean(w))**2) ) #equation (2)

• Right, I have only considered $(X^\top X)^{-1}$ but not $x_v^\top(X^{\top}X)^{-1}x_v$, or? – schn Feb 1 at 19:40
• Also, I guess $(X^\top X)^{-1}$ is not a matrix in my example, since $x_v=(x_0^1-\bar{x}_{1\cdot})$. – schn Feb 1 at 19:46
• It's a bit strange, since $(X^\top X)^{-1}= \begin{pmatrix} \frac{1}{n} & 0 \\ 0 & \frac{1}{\sum_{i=1}^n (x_i-\bar{x})^2 }\end{pmatrix}$ – schn Feb 1 at 20:09
• schn, That's unlikely and it's certainly not the case for the x in your question. See for yourself: compute it as solve(crossprod(cbind(1,x))) (or use w in place of x if that's your intent). – whuber Feb 2 at 12:55