sp <- gsub("setosa", 0, iris$Species)
sp <- gsub("versicolor", 1, sp)
iris$Species <- as.numeric(gsub("virginica", 2, sp))
iris$var <- rep(c(0,1,2), nrow(iris)/3)
fit <- lm (Sepal.Length ~ factor(Species):factor(var)-1, data = iris)

Call:
lm(formula = Sepal.Length ~ factor(Species):factor(var) - 1,
    data = iris)

Residuals:
     Min       1Q   Median       3Q      Max
-1.54706 -0.31875 -0.05165  0.31618  1.32941

Coefficients:
                              Estimate Std. Error t value Pr(>|t|)
factor(Species)0:factor(var)0   5.0529     0.1249   40.47   <2e-16 ***
factor(Species)1:factor(var)0   5.7706     0.1249   46.22   <2e-16 ***
factor(Species)2:factor(var)0   6.7563     0.1287   52.50   <2e-16 ***
factor(Species)0:factor(var)1   5.0118     0.1249   40.14   <2e-16 ***
factor(Species)1:factor(var)1   6.0187     0.1287   46.77   <2e-16 ***
factor(Species)2:factor(var)1   6.4471     0.1249   51.64   <2e-16 ***
factor(Species)0:factor(var)2   4.9500     0.1287   38.46   <2e-16 ***
factor(Species)1:factor(var)2   6.0235     0.1249   48.24   <2e-16 ***
factor(Species)2:factor(var)2   6.5706     0.1249   52.62   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5148 on 141 degrees of freedom
Multiple R-squared:  0.9928,    Adjusted R-squared:  0.9924
F-statistic:  2174 on 9 and 141 DF,  p-value: < 2.2e-16

lsmeans(fit, pairwise ~ factor1 | factor2)