data <- scan(what = list(pie="",hours=0,taste=0), multi.line=F)
 custard        4.1       5.5
 custard        3.2       4.9
 custard        2.3       4.6
 custard        3.0       5.1
 custard        2.9       4.8
 custard        3.0       4.8
 custard        4.1       5.3
 custard        2.1       4.6
 custard        2.5       4.6
 custard        4.2       5.5
 custard        3.1       5.0
 custard        4.2       5.6
 mincemeat      3.7       5.6
 mincemeat      4.5       6.2
 mincemeat      4.2       6.2
 mincemeat      5.4       7.4
 mincemeat      4.4       6.4
 mincemeat      3.5       5.7
 mincemeat      5.0       6.8
 mincemeat      4.9       6.7
 mincemeat      3.0       5.1
 mincemeat      3.8       5.9
 mincemeat      3.3       5.2
 mincemeat      3.0       5.1
 pumpkin        5.9       6.4
 pumpkin        4.3       5.8
 pumpkin        4.7       5.6
 pumpkin        6.3       6.5
 pumpkin        6.4       6.7
 pumpkin        5.1       5.9
 pumpkin        5.6       6.3
 pumpkin        5.2       5.8
 pumpkin        4.3       5.6
 pumpkin        4.0       5.3
 pumpkin        4.4       5.6
 pumpkin        4.1       5.7

data 

dessert <- data.frame(data)
rm(data)

dessert$pie  <- as.factor(dessert$pie)

 # look at the data

plot(dessert$hours,dessert$taste,pch=as.numeric(dessert$pie),
      col=as.numeric(dessert$pie),ylab="Taste",xlab="Hours")
legend("bottomright",levels(dessert$pie),pch=1:3,col=1:3)

abline(lm(taste[pie=="custard"] ~ hours[pie=="custard"], data=dessert),lwd=2,col="black")
abline(lm(taste[pie=="mincemeat"] ~ hours[pie=="mincemeat"], data=dessert),lwd=2,col="red")
abline(lm(taste[pie=="pumpkin"] ~ hours[pie=="pumpkin"], data=dessert),lwd=2,col="green")


 # check equal slopes assumption

dessert.lm1 <- lm(taste ~ hours +pie +hours:pie, data=dessert)
anova(dessert.lm1)
summary(dessert.lm1)


 # find median and quartiles of the covariate

summary(dessert$hours)

 # calculate adjusted means at the median and quartiles of the covariate

library(effects)

effect("pie",dessert.lm1,given.values=c(hours=c(3.2)))
effect("pie",dessert.lm1,given.values=c(hours=c(3.2)))$se

effect("pie",dessert.lm1,given.values=c(hours=c(4.2)))
effect("pie",dessert.lm1,given.values=c(hours=c(4.2)))$se

effect("pie",dessert.lm1,given.values=c(hours=c(4.8)))
effect("pie",dessert.lm1,given.values=c(hours=c(4.8)))$se


 # pairwise multiple comparison tests at different x values,
 # using the multcomp package

library(multcomp)

 # at x = 3.2

contrdiff3.2 <- rbind("cust -mince at 3.2"  <- c(0,0,-1,0,-3.2,0),
                    "mince - pump at 3.2"     <- c(0,0,1,-1,3.2,-3.2),
                    "cust - pump at 3.2"      <- c(0,0,0,-1,0,-3.2))
   
summary(glht(dessert.lm1, linfct= contrdiff3.2))

 # at x = 4.2

contrdiff4.2 <- rbind("cust -mince at 4.2"  <- c(0,0,-1,0,-4.2,0),
                   "mince - pump at 4.2"     <- c(0,0,1,-1,4.2,-4.2),
                   "cust - pump at 4.2"      <- c(0,0,0,-1,0,-4.2))
 
summary(glht(dessert.lm1, linfct= contrdiff4.2))


# at x = 4.8:

contrdiff4.8 <- rbind("cust -mince at 4.8"  <- c(0,0,-1,0,-4.8,0),
                   "mince - pump at 4.8"     <- c(0,0,1,-1,4.8,-4.8),
                   "cust - pump at 4.8"      <- c(0,0,0,-1,0,-4.8))
 
summary(glht(dessert.lm1, linfct= contrdiff4.8))