1
$\begingroup$

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?

$\endgroup$
2
  • $\begingroup$ 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$. $\endgroup$
    – Fiodor1234
    Feb 1, 2021 at 19:17
  • $\begingroup$ 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. $\endgroup$
    – psie
    Feb 1, 2021 at 19:33

1 Answer 1

1
$\begingroup$

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)
$\endgroup$
4
  • $\begingroup$ Right, I have only considered $(X^\top X)^{-1}$ but not $x_v^\top(X^{\top}X)^{-1}x_v$, or? $\endgroup$
    – psie
    Feb 1, 2021 at 19:40
  • $\begingroup$ Also, I guess $(X^\top X)^{-1}$ is not a matrix in my example, since $x_v=(x_0^1-\bar{x}_{1\cdot})$. $\endgroup$
    – psie
    Feb 1, 2021 at 19:46
  • $\begingroup$ 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} $ $\endgroup$
    – psie
    Feb 1, 2021 at 20:09
  • $\begingroup$ 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). $\endgroup$
    – whuber
    Feb 2, 2021 at 12:55

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.