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, Parallel Lines Model")

parallel.lines.model <- lm(flowers~intensity + time)

## to see the regression output, go to the command line, type:
## summary(parallel.lines.model)

## to see the coefficients alone, go to the command line, type:
## parallel.lines.model$coefficients
## We will refer to the components of this Splus object to
## make the plot.

## fitted coefficients from parallel lines model
beta.0 <- parallel.lines.model$coefficients[1]
beta.intensity <- parallel.lines.model$coefficients[2]
beta.time <- parallel.lines.model$coefficients[3]

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

## add a line for the timing=1 group
abline((beta.0+beta.time), beta.intensity, lty=2)
## "lty=2" denotes a dotted line
## default is "lty=1" a solid line