#example with two categories; one indicator variable
fert.data2=read.table("FertVariety2.txt",header=TRUE)
attach(fert.data2)

n=length(yield)
#define variables for the model
#response:
y=yield
#explanatory variables:
x1=fert
#indicator for variety
#start with a vector of n zeroes
x2=rep(0,n)
#if variety is A, then x2=1
x2[variety=="A"]=1
#and the interaction is simply the product:
x3=x1*x2
#fit the model:
lm(y~x1+x2+x3)

#example with three categories; two indicator variable
fert.data3=read.table("FertVariety3.txt",header=TRUE)
attach(fert.data3)
#define variables for the model
#compute sample size
n=length(yield)
#response:
y=yield
#explanatory variables:
x1=fert
#indicator for variety, will need two variables
#start with a vector of n zeroes
x2=rep(0,n)
#if variety is A, then x2=1
x2[variety=="A"]=1
#repeat for variety B
x3=rep(0,n)
#if variety is A, then x2=1
x3[variety=="B"]=1

#interaction given by products:
x4=x1*x2
x5=x1*x3

#fit the model:
lm(y~x1+x2+x3+x4+x5)