Getting data

y <- c(91,105,106,108,88,91,58,82,81,65,61,48,61,43,33,36)
x1 <- c(16.7,17.1,18.2,18.1,17.2,18.2,16,17.2,18,17.2,16.9,17.1,18.2,17.3,17.5,16.6)
x2 <- c(1:16)

Getting least square of including two factor interaction with using lm function

model1 <- lm(y~x1+x2+x1*x2)
model1
## 
## Call:
## lm(formula = y ~ x1 + x2 + x1 * x2)
## 
## Coefficients:
## (Intercept)           x1           x2        x1:x2  
##   -220.2707      19.1985      11.1880      -0.9206

Getting least square estimate without lm .

one <- c(rep(1,16))
X <- cbind(one,x1,x2,x1*x2)
Y <- cbind(y)
firstTerm <- ginv(t(X) %*% X)
secondTerm <- t(X) %*% Y
finalAns <- firstTerm %*% secondTerm
finalAns
##                y
## [1,]  -0.7854925
## [2,]   6.5328127
## [3,] -11.3158845
## [4,]   0.3808286