2
d = data.frame(
    Temperature = c(rep("Cool", 6), rep("Warm", 6)),
    Bact = c(rep("Bact 1", 2), rep("Bact 2", 2), rep("Bact 3", 2), rep("Bact 1", 2), rep("Bact 2", 2), rep("Bact 3", 2)),
    Time = c(15.23,14.32,14.77,15.12,14.05,15.48,14.13,16.13,16.44,14.82,17.96,16.65)
)

I self-created a small data frame for a two-way ANOVA. I want to perform a two-way ANOVA model by

summary(aov(Time~Bact*Temperature, data=d))

Time is the dependent variable while Bact and Temperature are two categorical independent variables.

Instead of doing it in the ANOVA way, I want to learn and prove ANOVA can also be done with a linear regression model. I want to convert my data into dummy variables and perform a linear regression on it. I expect I'll recover the same results. The dummy variables will also include interactions effects between Bact and Temperature.

The problem is that, I don't know how to convert my data frame into dummy variables such that it can be used in the lm() function.

LyzandeR
  • 37,047
  • 12
  • 77
  • 87
ABCD
  • 7,914
  • 9
  • 54
  • 90

2 Answers2

2

I kind of do the same with you. I want to feel in control so whenever I have time I design the dummies myself with the following:

d = data.frame(
  Temperature = c(rep("Cool", 6), rep("Warm", 6)),
  Bact = c(rep("Bact 1", 2), rep("Bact 2", 2), rep("Bact 3", 2), rep("Bact 1", 2), rep("Bact 2", 2), rep("Bact 3", 2)),
  Time = c(15.23,14.32,14.77,15.12,14.05,15.48,14.13,16.13,16.44,14.82,17.96,16.65)
)

which is:

> d
   Temperature   Bact  Time
1         Cool Bact 1 15.23
2         Cool Bact 1 14.32
3         Cool Bact 2 14.77
4         Cool Bact 2 15.12
5         Cool Bact 3 14.05
6         Cool Bact 3 15.48
7         Warm Bact 1 14.13
8         Warm Bact 1 16.13
9         Warm Bact 2 16.44
10        Warm Bact 2 14.82
11        Warm Bact 3 17.96
12        Warm Bact 3 16.65

So you only need to dummify factors (temperature,bact) so the following process works:

xfactors <- Filter(is.factor,d) #filter only the factors to dummify
b <- data.frame(matrix(NA,nrow=nrow(xfactors),ncol=1)) #make empty data.frame to initiate b
for ( i in 1:ncol(xfactors)) { #start loop
  a <- data.frame(model.matrix(~xfactors[,i])) #make dummies here
  b <- cbind(b, a[-1]) #remove intercept and combine dummies
}
b <- data.frame(b[-1]) #make a data.frame
#the reference dummy gets excluded automatically by model.matrix
colnames(b) <- c('warm' , 'bact2' , 'bact3') #you will probably want to change the names to sth smaller

> b
   warm bact2 bact3
1     0     0     0
2     0     0     0
3     0     1     0
4     0     1     0
5     0     0     1
6     0     0     1
7     1     0     0
8     1     0     0
9     1     1     0
10    1     1     0
11    1     0     1
12    1     0     1

Then to run the model:

new_data <- cbind(b, Time=d$Time) #add time to the data
mymod <- lm(Time ~ warm*bact2+warm*bact3, data=new_data) #compute lm with interactions
#you shouldn't compute the interactions between dummy variables because they come from the same variable

Output:

> summary(mymod)

Call:
lm(formula = Time ~ warm * bact2 + warm * bact3, data = new_data)

Residuals:
   Min     1Q Median     3Q    Max 
 -1.00  -0.67   0.00   0.67   1.00 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  14.7750     0.6873  21.498 6.61e-07 ***
warm          0.3550     0.9719   0.365    0.727    
bact2         0.1700     0.9719   0.175    0.867    
bact3        -0.0100     0.9719  -0.010    0.992    
warm:bact2    0.3300     1.3745   0.240    0.818    
warm:bact3    2.1850     1.3745   1.590    0.163    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9719 on 6 degrees of freedom
Multiple R-squared:  0.6264,    Adjusted R-squared:  0.3151 
F-statistic: 2.012 on 5 and 6 DF,  p-value: 0.2097
LyzandeR
  • 37,047
  • 12
  • 77
  • 87
  • 1
    `warm * bact2` == `warm + bact2 + warm:bact2`. Adding another `warm*bact3` will indeed add `warm` twice but R is smart enough to exclude multiple instances of variables with the same name and include them only once. Try `lm(Time ~ warm + warm + warm, data=new_data)` and see how many coefficients you ll get :) – LyzandeR Jan 09 '15 at 12:03
0

lm() will create dummy variables for you. No need to create them yourself:

m <- lm(Time ~ Bact*Temperature, data = d)
anova(m)

Edit

If you want to peer under the hood of lm(), you can see the design matrix with model.matrix(m)

davechilders
  • 8,693
  • 2
  • 18
  • 18
  • This is a self-study problem. I'm interested in learning the statistics and doing the dirty work myself otherwise hidden by the R functions. – ABCD Jan 09 '15 at 05:22
  • @StudentT If you want to generate the "X" matrix, you could use `model.matrix(Time ~ Bact*Temperature, data = d)` but then you can't use such a matrix as an input to `lm` anymore. It's weird to ask R just to do half the work for you. – MrFlick Jan 09 '15 at 06:20