Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I have two datasets, A and B, of weighted (x,y) pairs. I computed the best fit lines, L_A and L_B, respectively, of these datasets, and then computed the intersection of these two lines, (x*,y*).

Now, A and B are drawn from random distributions (of the form y=mx+b+e, where e is a normally distributed error term). I'd like to compute a confidence interval around x*. How can I do this?

A naive approach I tried involved computing confidence intervals for the slope/intercept estimates of L_A and L_B, generating approximate distributions for those 4 parameters, and repeatedly randomly sampling from these distributions. But a given line's slope and intercept estimates are not independent of each other, and I wasn't sure how to salvage this approach to account for that.

share|improve this question
add comment

1 Answer

up vote 7 down vote accepted

One straightforward way is to obtain the maximum likelihood estimator of $x^{*}$ directly. Using the first subscript to designate the line ($1$ or $2$), the model is

$$y_{ij} = m_i(x_{ij} - x^{*}) + y^{*} + \varepsilon_{ij},$$

$1 \le j \le n_i$, $\varepsilon_{ij} \sim \text{Normal}(0, \sigma_i / \omega_{ij})$ and independent. The $n_i$ count the data for each line. The parameters are the point of intersection $(x^{*}, y^{*})$, the slopes $m_1$ and $m_2$, and (as nuisance parameters) the scale factors $\sigma_1$ and $\sigma_2$. The $\omega_{ij}$ are specified weights (not parameters). The standard ML machinery will provide a confidence interval for any of the parameters, including $x^{*}$.

Two lines with confidence bands

This illustration show two linear fits from $10$ independent samples of two lines, each passing through $(7,10)$, with slopes $1$ and $1/3$. The colored bands are $95\%$ confidence bands around the fits (using least squares). The estimated point of intersection, shown as a large black dot, occurs at $(6.3, 9.4)$. A $95\%$ confidence interval for its x-coordinate is portrayed as a dashed black segment: it extends from $4.4$ to $8.3$.


Some details

(added in response to comments)

This model is identical to performing two separate weighted least squares regressions. It merely combines their 3+3 = 6 parameters (intercept, slope, and scale) in a way that isolates the x-coordinate of the point of intersection. Therefore the parameter estimates will be the same. The point is that in fitting the combined model, any maximum likelihood procedure will report standard errors (equivalently, confidence intervals) for the parameters and this is how we solve the original problem.

I would also like to point out that the there is some freedom in choosing the weights (which can be helpful in persuading an optimization routine to behave well). The foregoing analysis shows it suffices to examine the case of fitting a single line. Writing $\omega_i$ for the weights the negative log likelihood is equal to

$$-\log(\Lambda) = \frac{1}{2}\sum_i{\log(2 \pi (\sigma/\omega_i)^2) + \frac{(m(x_i-x^*) + y^* - y_i)^2}{(\sigma/\omega_i)^2}}.$$

Removing additive constants (which don't affect the ML procedure) this simplifies somewhat to

$$\sum_i{\log(\sigma) + \frac{u_i^2\omega_i^2}{2\sigma^2}}$$

where $u_i^2 = (m(x_i-x^*) + y^* - y_i)^2$. To estimate $\sigma$ we equate the derivative with $0$:

$$0 = \frac{\partial(-\log(\Lambda))}{\partial \sigma} = \sum_i^n{\frac{1}{\sigma} - \frac{u_i^2\omega_i^2}{\sigma^3}},$$

implying

$$(\hat{\sigma})^2 = \frac{1}{n}\sum_i^n{u_i^2 \omega_i^2}.$$

This shows that the estimate of $\sigma^2$ is directly proportional to the scale of the weights. This is clear: if we multiply all weights by a positive value $\lambda$, then we get exactly the same model by dividing $\sigma$ by $\lambda$. Consequently, we are free to normalize the weights as we wish. A good choice is to make $\sum{\omega_i^2}=n$, because this is what we would get for an unweighted model. In particular, note that the actual value of $\sigma$ is meaningless except in comparison to the $L^2$ norm of the weights. Accordingly, the standard error of $\hat{\sigma}$ is also meaningless except in comparison to the weights.

share|improve this answer
1  
Could you expand upon how to use "standard ML machinery" to derive the desired confidence interval? Do I need to get my hands dirty with some algebra and calculus, or is there a simpler approach? –  David Sep 13 '11 at 19:34
    
@David You can do this with calculus and a modest understanding of optimization software, but the whole point is that there's loads of software around that will do everything for you. Search RSeek for "maximum likelihood" or "MLE" for instance. If you have access to commercial software like Stata or Mathematica, routines for computing ML confidence intervals are built in. –  whuber Sep 13 '11 at 21:30
    
Thank you, whuber, this was very helpful and educational. I derived the log-likelihood as a function of the parameters, and passed the function into R via nlm(). However, nlm() doesn't seem up to task, as its parameter estimates appear to vary greatly depending on the starting parameter values I pass in. Additionally, the estimates are fairly off from the intersection of the 2 best fit lines which I computed directly. I'm not sure if this is due to nlm()'s shortcomings or if I'm not using it properly. Any suggestions? –  David Sep 14 '11 at 17:43
    
@David nlm() does have shortcomings, but they shouldn't be that bad! There's a reliable way to give it good starting values: fit the lines separately using weighted OLS and compute their intersection. The values you get for $x^{*}, y^{*}, m_1, m_2$ should already be optimal for nlm(). It will go through one iteration to compute the gradient and hessian, using those to produce confidence intervals. That's really all you're using it for. If it gets different answers than you started with, then there's probably something wrong with your code. –  whuber Sep 14 '11 at 18:46
1  
Thank you so much! The two-pass method worked, resulting in exact matches between the nlm() values and the values obtained by computing the lines and intersection directly. The SE value determined from the hessian seems reasonable. I am very grateful for your help! –  David Sep 15 '11 at 18:35
show 7 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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