What if my linear regression data contains several co-mingled linear relationships? Let's say I am studying how daffodils respond to various soil conditions.  I have collected data on the pH of the soil versus the mature height of the daffodil.  I'm expecting a linear relationship, so I go about running a linear regression.
However, I didn't realize when I started my study that the population actually contains two varieties of daffodil, each of which responds very differently to soil pH.  So the graph contains two distinct linear relationships:

I can eyeball it and separate it manually, of course.  But I wonder if there is a more rigorous approach.
Questions:


*

*Is there a statistical test to determine whether a data set would be better fit by a single line or by N lines?

*How would I run a linear regression to fit the N lines?  In other words, how do I disentangle the co-mingled data?
I can think of some combinatorial approaches, but they seem computationally expensive.

Clarifications:


*

*The existence of two varieties was unknown at the time of data collection.  The variety of each daffodil was not observed, not noted, and not recorded.

*It is impossible to recover this information.  The daffodils have died since the time of data collection.
I have the impression that this problem is something similar to applying clustering algorithms, in that you almost need to know the number of clusters before you start.  I believe that with ANY data set, increasing the number of lines will decrease the total r.m.s. error.  In the extreme, you can divide your data set into arbitrary pairs and simply draw a line through each pair.  (E.g., if you had 1000 data points, you could divide them into 500 arbitrary pairs and draw a line through each pair.)  The fit would be exact and the r.m.s. error would be exactly zero.  But that's not what we want.  We want the "right" number of lines.
 A: I think Demetri's answer is a great one if we assume that you have the labels for the different varieties.  When I read your question that didn't seem to be the case to me.  We can use an approach based on the EM algorithm to basically fit the model that Demetri suggests but without knowing the labels for the variety. Luckily the mixtools package in R provides this functionality for us. Since your data is quite separated and you seem to have quite a bit it should be fairly successful.
library(mixtools)

# Generate some fake data that looks kind of like yours
n1 <- 150
ph1 = runif(n1, 5.1, 7.8)
y1 <- 41.55 + 5.185*ph1 + rnorm(n1, 0, .25)

n2 <- 150
ph2 <- runif(n2, 5.3, 8)
y2 <- 65.14 + 1.48148*ph2 + rnorm(n2, 0, 0.25)

# There are definitely better ways to do all of this but oh well
dat <- data.frame(ph = c(ph1, ph2), 
                  y = c(y1, y2), 
                  group = rep(c(1,2), times = c(n1, n2)))

# Looks about right
plot(dat$ph, dat$y)

# Fit the regression. One line for each component. This defaults
# to assuming there are two underlying groups/components in the data
out <- regmixEM(y = dat$y, x = dat$ph, addintercept = T)

We can examine the results
> summary(out)
summary of regmixEM object:
          comp 1    comp 2
lambda  0.497393  0.502607
sigma   0.248649  0.231388
beta1  64.655578 41.514342
beta2   1.557906  5.190076
loglik at estimate:  -182.4186 

So it fit two regressions and it estimated that 49.7% of the observations fell into the regression for component 1 and 50.2% fell into the regression for component 2. The way I simulated the data it was a 50-50 split so this is good.
The 'true' values I used for the simulation should give the lines:
y = 41.55 + 5.185*ph and y = 65.14 + 1.48148*ph
(which I estimated 'by hand' from your plot so that the data I create looks similar to yours) and the lines that the EM algorithm gave in this case were:
y = 41.514 + 5.19*ph and y = 64.655 + 1.55*ph
Pretty darn close to the actual values.
We can plot the fitted lines along with the data
plot(dat$ph, dat$y, xlab = "Soil Ph", ylab = "Flower Height (cm)")
abline(out$beta[,1], col = "blue") # plot the first fitted line
abline(out$beta[,2], col = "red") # plot the second fitted line


A: I'll focus on the question of statistical significance since Dason already covered the modeling part.
I am unfamiliar with any formal tests for this (which I am sure exist), so I'll just throw some ideas out there (and I'll probably add R code and technical details later).
First, it is convenient to infer the classes. Presuming you have two lines fit to the data, you can approximately reconstruct the two classes by assigning each point to the class of the line closest to it. For points near the intersection, you will run into issues, but for now just ignore those (there may be a way to get around this, but for now just hope that this won't change much).
The way to do this is to choose $x_{l}$ and $x_{r}$ (soil pH values) with $x_{l} \leq x_{r}$ such that the parts left of $x_{l}$ are sufficiently separated and the parts right of $x_{r}$ are sufficiently separated (the closest point where the distributions don't overlap).
Then there are two natural ways I see to go about doing this.
The less fun way is to just run your original dataset combined with the inferred class labels through a linear regression as in Demetri's answer.
A more interesting way to do so would be through a modified version of ANOVA.
The point is to create an artificial dataset that represents the two lines (with similar spread between them) and then apply ANOVA. Technically, you need to do this once for the left side, and once for the right (i.e. you'll have two artificial datasets).
We start with the left, and apply a simple averaging approach to get two groups. Basically, each point in say the first class is of the form
$$ y^{(i)}_{1} = \beta_{1,1} x_{1}^{(i)} + \beta_{1,0} + e_{1}^{(i)}$$
so we are going to replace the linear expression
$\beta_{1,1} x_{1}^{(i)} + \beta_{1,0}$
by a constant, namely the average value of the linear term or
 $$ \beta_{1,1} x^{\mathrm{avg}} + \beta_{1, 0}$$
where $x^{\mathrm{avg}}_{l}$ is literally the average $x$ value for the left side (importantly, this is over both classes, since that makes things more consistent). That is, we replace $y_{1}^{(i)}$ with
  $$ \tilde{y}_{1}^{(i)} = \beta_{1,1} x^{\mathrm{avg}} + \beta_{1, 0} + e_{1}^{(i)},$$
and we similarly for the second class. That is, your new dataset consists of the collection of $\tilde{y}_{1}^{(i)}$ and similarly $\tilde{y}_{2}^{(i)}$.
Note that both approaches naturally generalize to $N$ classes.
A: Is it possible that including both in the same chart is an error?  Given that the varieties behave completely different is there any value in overlapping the data?  It seems to me that you are looking for impacts to a species of daffodil, not the impacts of similar environments on different daffodils.  If you have lost the data that helps determine species "A" from species "B" you can simply group behavior "A" and behavior "B" and include the discovery of two species in your narrative.  Or, if you really want one chart, simply use two data sets on the same axis. I don't have anywhere near the expertise that I see in the other responses given so I have to find less "skilled" methods. I would run a data analysis in a worksheet environment where the equations are easier to develop.  Then, once the groupings become obvious, create the two separate data tables followed by converting them into charts/graphs.  I work with a great deal of data and I often find that my assumptions of differing correlations turn out wrong; that is what data is supposed to help us discover.  Once I learn that my assumptions are wrong, I display the data based upon the behaviors discovered and discuss those behaviors and resulting statistical analyses as part of the narrative.  
