plot(intensity,flowers,type="n",xlab="Intensity",ylab="Flowers")
points(intensity[time==0],flowers[time==0],pch="0")
points(intensity[time==1],flowers[time==1],pch="1")
title("Flowers vs. Intensity, Separate Lines Model")

## fit the separate lines model
separate.lines.model <- lm(flowers~intensity * time)

## note this can also be fit with:
## separate.lines.model <- lm(flowers~intensity + time + intensity:time)

beta.0 <- separate.lines.model$coefficients[1]
b.intensity <- separate.lines.model$coefficients[2]
beta.time <- separate.lines.model$coefficients[3]
beta.intensity.time <- separate.lines.model$coefficients[4]

## add the regression line for timing=0
abline(beta.0,beta.intensity)

## add the regression line for timing=1
abline((beta.0 + beta.time),(beta.intensity + beta.intensity.time), lty=2)