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

Question

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?

Answer 1

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)