1
$\begingroup$

A simple question I hope.

I have an experimental design where I measure some response (let's say blood pressure) from two groups: a control group and an affected group, where both are given three treatments: t1, t2, t3. The data are not paired in any sense.

Here is an example data:

set.seed(1)
df <- data.frame(response = c(rnorm(5,10,1),rnorm(5,10,1),rnorm(5,10,1),
                         rnorm(5,7,1),rnorm(5,5,1),rnorm(5,10,1)),
                 group = as.factor(c(rep("control",15),rep("affected",15))),
                 treatment = as.factor(rep(c(rep("t1",5),rep("t2",5),rep("t3",5)),2)))

What I am interested in is quantifying the effect that each treatment has on the affected group relative to the control group. How would I model this, say using an linear model (for example lm in R)?

Am I wrong thinking that:

lm(response ~ 0 + treatment * group, data = df)

which is equivalent to:

lm(response ~ 0 + treatment + group + treatment:group, data = df)

is not what I need? I think that in this model the treatment:group interaction terms are relative to the mean over all baseline group and baseline treatment measurements.

I therefore thought that this model:

lm(response ~ 0 + treatment:group, data = df)

is what I need but it's quantifying each combination of treatment and group interaction terms: treatmentt1:groupcontrol treatmentt1:groupaffected treatmentt2:groupcontrol treatmentt2:groupaffected treatmentt3:groupcontrol treatmentt3:groupaffected

So perhaps this model:

lm(response ~ 0 + treatment + treatment:group, data = df)

is the correct one?

Although in addition to quantifying each combination of treatment and groupaffected interaction term it's also quantifying the effect of each treatment. I'm not sure what is the baseline each of the treatment and groupaffected interaction terms are compared to in this model.

Help would be appreciated.

Also, let's say I ran a fourth treatment which is actually the combination of two treatments, say t1+t3, where I don't know what the expectation of their combined effect is: additive/subtractive or synergistic. Is there any way this can be combined?

$\endgroup$
2
$\begingroup$

You should not exclude the intercept. Instead, you should make sure that the control is the reference group for treatment contrasts:

df$group <- relevel(df$group, "control")

Then you should of course plot your data:

library(ggplot2)
ggplot(df, aes(x = treatment, y = response, color = group)) +
  geom_boxplot()

resulting plot

It looks like the control group doesn't differ between treatments and the affected group seems to exhibit lower values for treatments t1 and t2, but not t3.

I would now fit the following model:

fit <- lm(response ~ group * treatment, data = df)
summary(fit)
#Call:
#lm(formula = response ~ group * treatment, data = df)
#
#Residuals:
#    Min      1Q  Median      3Q     Max 
#-2.2528 -0.4973  0.1965  0.5250  1.4737 
#
#Coefficients:
#                           Estimate Std. Error t value Pr(>|t|)    
#(Intercept)               10.129270   0.438666  23.091  < 2e-16 ***
#groupaffected             -2.669703   0.620367  -4.303 0.000244 ***
#treatmentt2                0.005866   0.620367   0.009 0.992534    
#treatmentt3               -0.091147   0.620367  -0.147 0.884419    
#groupaffected:treatmentt2 -2.384202   0.877331  -2.718 0.012013 *  
#groupaffected:treatmentt3  2.283003   0.877331   2.602 0.015626 *  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 0.9809 on 24 degrees of freedom
#Multiple R-squared:  0.8226,   Adjusted R-squared:  0.7857 
#F-statistic: 22.26 on 5 and 24 DF,  p-value: 2.691e-08

Here, the intercept represents the mean of the control group at t1. We see that the control groups at t2 and t3 don't differ significantly from this mean (as the treatment coefficients tell us). We also see that we have a significant difference between the groups at t1 and an even larger difference at t2, but approximately no difference (-2.67 + 2.28) at t3. You could test the significance of the latter using a post-hoc test, but it is probably sufficient to look at the least squares means:

library(lsmeans)
lsmeans(fit, c("group", "treatment"))
#group    treatment    lsmean        SE df lower.CL  upper.CL
#control  t1        10.129270 0.4386656 24 9.223909 11.034631
#affected t1         7.459567 0.4386656 24 6.554206  8.364928
#control  t2        10.135136 0.4386656 24 9.229774 11.040497
#affected t2         5.081231 0.4386656 24 4.175869  5.986592
#control  t3        10.038123 0.4386656 24 9.132762 10.943484
#affected t3         9.651423 0.4386656 24 8.746062 10.556784
#
#Confidence level used: 0.95 

As you see, the confidence intervals of the groups at t3 overlap strongly.

Obviously, you need to employ the usual diagnostics to check for variance heterogeneity, influential points etc.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ If don't exclude the intercept in lm(response ~ 0 + treatment + treatment:group, data = df), meaning: lm(response ~ treatment + treatment:group, data = df) I get all 3 interaction terms: treatmentt1:groupaffected, treatmentt2:groupaffected, and treatmentt3:groupaffected quantifying the effect of each treatment on the affected group relative to the control group. That seems to me more intuitive than the summary of: lm(response ~ group * treatment, data = df) $\endgroup$ – dan Nov 12 '15 at 17:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.