Fitting data based on an unobserved variable

Question

I have pairs (x, y) which I would like to regress (x independent, y dependent). Plotting them, I see distinct bands which can be attributed to a third variable. Unfortunately this variable cannot be directly measured.

Is there a good way to both fit the data to an equation and determine the value of the binary variable (for each x)? The equation I was considering is

$$y=ax^b(1+cz)$$

where $z$ is the unobserved binary variable. (Of course there are two ways to fit the data depending on which group is chosen as 0 and which as 1, but restricting $c>0$ should make the choice essentially unique.)

Unless there's reason to prefer something else I would look for a least-squares fit.

Answer 1

Here is a solution that is close to least squares, using a Gaussian mixture model. We use a family of probability distributions with PDFs that capture both centers of the conditional distribution of $y|x$:

$$p(y,x; (a,b,\delta,\pi, \sigma)) =(1-\pi)\phi_\sigma(y-a x^b)+\pi\phi_\sigma(y - a x^b(1+\delta)).$$

In this notation $\phi_\sigma$ is the density function for a Normal$(0, \sigma)$ distribution. The weighted linear combination captures the idea of two vertically separated bands (as a mixture, with the chance of being in the upper band being $\pi$). There is Normal variation around the centers at $a x^b$ and $a x^b(1+\delta)$.

This family depends on three parameters of interest, $a$ and $b$ (describing the functional form) and $\delta$ (expressing the separation between the two bands). It uses two ancillary parameters, $\pi$ (the proportions in each band, assumed constant across $x$) and $\sigma$ (the residual standard deviation).

To illustrate, here are plots of $p(y, 1/6; (1/2, 1, 1, 1/3, 1/3))$ (blue, to the left) and $p(y, 1/2; (2, 1/2, 1, 1/3, 1/3))$ (red, to the right):

Fitting these with maximum likelihood is the best analog of least squares available. As an example, consider these data generated by evaluating $a x_i^b(1 + \delta u_i)+\varepsilon_i$ for $a=2$, $b=1/2$, and $\delta=1$ at $x_i = 1/64, 2/64, \ldots, 1$ where $u_i$ were drawn independently from a Bernoulli$(1/3)$ distribution and, independently of those, $\varepsilon_i$ were independently drawn from a Normal$(0, \sigma)$ distribution, $\sigma = 1/3$:

This is not a lot of data and the separation into two bands is not completely clear, so it forms a moderately severe test of the method. The MLE gives

$$\hat{a} = 1.99, \hat{b} = 0.509, \hat{\delta} = 1.07, \hat{\pi} = 0.27, \hat{\sigma} = 0.33$$

compared to

$$a = 2.00, b = 0.50, \delta = 1.00, \pi = 0.33, \sigma = 0.33.$$

Apart from $\pi$, which will be difficult to estimate with such few data and is of little consequence, these estimates are extremely close to the actual values used to generate the data, and this is no accident: repeated simulations continue to get good estimates. These fits agree nicely with the data:

Answer 2

I'm not sure that you conceptually can fit a parametrized curve directly on unobserved data. Maybe a two-step approach would help. but I admit, I have never seen or done such a thing: First fit the curve only for 2 variables $x$ and $y$. If the bands of points you mentioned do not cross then the fitted curve should lie somewhere in between the two bands. Then use the residuals to "reconstruct" the third variable $z$ by setting it to 0 if the residual is negative, say, and setting it to 1 if the residual is positive. Then use the variable $z$ together with $x$ and $y$ in your proposed equation. In the case of a linear regression this R code could demonstrate this procedure:

# x-values
x1<-seq(0,10,by=0.1)

# the two groups of values lying around two lines
y1<-2+3*x1+rnorm(101,0,1)
y2<-6+3*x1+rnorm(101,0,1)

# plot all points
x<-c(x1,x1)
y<-c(y1,y2)

plot(x,y)

# one model for both groups of points
lm.forboth<-lm(y~x)

# plot the fitted values
lines(x,lm.forboth$fitted.values)

# find those points whose residuals are >0 and let z be zero for them
# this is the reconstruction of z
posresid<-lm.forboth$residuals>0
z<-as.numeric(posresid)

# now use z to fit two lines
lm.separate<-lm(y~x+z)

# plot the fit of the two lines obtained by z=0,1
lines(x[posresid],lm.separate$fitted.values[posresid])
lines(x[!(posresid)],lm.separate$fitted.values[!(posresid)])

I hope this helps...

Answer 3

If you take logs, you would have

$$\log(y) = \log(a) + b\log(x) + \log(1+cz) + \epsilon.$$

Consider $\log(1+cz) + \epsilon$ as a random perturbation: call it $\eta$. If you regress $\log(y)$ on a column of "ones" and $\log(x)$, you get estimates of $\log(a)$, $b$ and the residuals $\hat\eta = \log(y) - \widehat\log(a) -\hat{b}x$.

If $\epsilon$ could be thought to have, for instance, a normal distribution, $\hat\eta$ would be roughly distributed as a mixture of normals with means $\log(1+c)$ apart. Then, you could use the EM algorithm as in here to obtain an estimate of the two means and their difference, along with the probabilities that $z$ takes value 0 or 1 conditional on the observed value of each residual.

Not given much thought, not tested: might be complete nonsense, please excuse me if it is.

Answer 4

This looks like a 'mixture of regressions' problem, where z is an unobserved random indicator variable denoting which curves each x-y observation belongs to. You would then, as F. Tusell suggests, use an EM algorithm to maximise the likelihood of all the regression parameters.

Here, EM would work by alternating the computation of the expected posterior probability of z for each x-y pair and then fitting the coefficients of the individual regressions on the assumption that these expectations were in fact the observed values of z. The regmixEM function from the R package mixtools will do this out of the box. There are other options if linear regression is not a good base model class.

The remaining task is to specify the individual regressions right. For this, you might fit the whole thing in log space as F. Tusell suggests, or add quadratic terms, or do whatever else makes sense in the context of the problem. This task is not obviously harder (or easier) than if z were observed and explicitly conditioned on.