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


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

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 ) 

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)

With output:

1.2313777318582835 1.8006635458783546
3.4259463217029915 6.053663015157194

In R and using predict.lm():




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?

  • $\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 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$ – schn Feb 1 at 19:33

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)
  • $\begingroup$ Right, I have only considered $(X^\top X)^{-1}$ but not $x_v^\top(X^{\top}X)^{-1}x_v$, or? $\endgroup$ – schn Feb 1 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$ – schn Feb 1 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$ – schn Feb 1 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 at 12:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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