# interpreting ancova and extracting estimated means

Imagine the following data:

ds <- data.frame(x=1:10,y=1:10,z=rep(c("A","B"),each=5))


the means for groups in z are:

library(plyr)
ddply(ds, "z", function(x) mean(x\$y))
#  z V1
#1 A  3
#2 B  8


do a couple of models:

m1 <- glm(y ~ x , data = ds)
m2 <- glm(y ~ z , data = ds)


In m1 the intercept is zero and the estimate for x is 1 as we would expect form the created data. Intercept meaning the value of y for x=0.

In m2 the intercept is 3 and the estimate for B is 5. This is because of the contrasts between the intercept with the estimate i.e. the estimate for B is 8. These results make sense in terms of the means of my data as calculated separately above.

so for m1 I interpret intercept as y where x=0 and slope and for m2 as this is an anova, intercept as mean of first level and the rest are differences between that level and the first.

However in the following models m3 and m4

m3 <- glm(y ~ x + z , data = ds)
m4 <- glm(y ~ x * z , data = ds)

# m3
#Coefficients:
#              Estimate Std. Error    t value Pr(>|t|)
#(Intercept)  1.123e-15  6.201e-16  1.812e+00    0.113
#x            1.000e+00  1.720e-16  5.814e+15   <2e-16 ***
#zB          -7.628e-17  9.880e-16 -7.700e-02    0.941

# m4
#Coefficients:
#              Estimate Std. Error    t value Pr(>|t|)
#(Intercept)  1.123e-15  8.714e-16  1.289e+00    0.245
#x            1.000e+00  2.627e-16  3.806e+15   <2e-16 ***
#zB          -9.349e-17  2.305e-15 -4.100e-02    0.969
#x:zB         3.130e-18  3.716e-16  8.000e-03    0.994


The coefficients for intercept and x are the same as m1 but how do I interpret the estimate for B and A as the estimate for B is -7.628e-17.

My question is how do I interpret the categorical estimates in m3 and from the output extract the estimated means, which I would report? and for m4 is the slope for x:zB 1.000e+00 + 3.130e-18 and for A 1.000e+00 - 3.130e-18?

Your models have a problem, you have added no error to your y so you get a perfect fit with your x variable. So let's fix this first:

set.seed(123)
ds <- data.frame(x=1:10,y=1:10+rnorm(10, 0, 0.5),z=rep(c("A","B"),each=5))
m1 <- glm(y ~ x , data = ds)
m2 <- glm(y ~ z , data = ds)
m3 <- glm(y ~ x + z , data = ds)
m4 <- glm(y ~ x * z , data = ds)
summary(m1)
summary(m2)
summary(m3)
summary(m4)


Now you get:

#m3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.3810     0.4124   0.924    0.386
x             0.9053     0.1144   7.914 9.77e-05 ***
zB            0.3547     0.6571   0.540    0.606


And:

#m4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.1552     0.4677  -0.332 0.751219
x             1.0840     0.1410   7.687 0.000254 ***
zB            2.3208     1.2374   1.876 0.109830
x:zB         -0.3575     0.1994  -1.793 0.123219


Now why do you get those odd estimates for B and for the interaction? If you look at R-squared in models m1 and m2 you'll notice that R^2 for m1 is 0.98 and for m2 is 0.77. That is, you are already explaining almost all the variability with just one variable so when you add another there's not much to add to the fit of the model. Look how the B coefficient is not significant in m3 but it was in m2. You can interpret the coefficient of B in m3 as the independent effect of B in y taking into account the effect of x in y. Since y is almost fully explained by x the estimated independent effect of B is much lower than what you'd expect.

• Thank you @Aghila +1 if I had more rep, I think that makes a lot of sense. However, if you could tell If I were to report the 2 means for z and 2 slopes forz in your m4 what they would be? So I know I am subrtactingadiing the effect sizes in the correct way. Dec 9, 2013 at 10:38