I want to analyse the effect of different treatment types (
control, treatment1, ..., treatment4) on the surface of specimens made of certain materials (
plastic, metal). The undamaged area of the surface is measured
after the treatment.
According to this design I specified a mixed model using lme4 as follows:
require("lme4") data <- read.csv("http://pastebin.com/raw/G4D8dh1f") mm1 <- lmer(undamaged_area ~ time*material*treatment + (1|specimen_id), data)
Is the mixed model the optimal choice in this case? I found some hints that an ANCOVA (something like
lm(undamaged_area_after ~ material*treatment + undamaged_area_before, data)) might be an alternative approach.
A closer look on the diagnostic plots of the mixed model makes me very suspicious:
plot(mm1); require("lattice"); qqmath(mm1)
Does the plots actually indicate a violation of model assumptions? Does the strange pattern come from a misspecification of the model?
Progress after donlelek's answer (=mixed model not required):
Just to be clear: The treatments were measured in different pieces of metal/plastic. So every piece is exactly measured twice - before and after the treatment. Thus, we are aiming for the ANOVA on damage, I guess. I had the impression to lose informations by just substracting the pre-post values. I did further research in the literature (with my limited knowledge in statistics). But according to "Pretest-posttest designs and measurement of change" the use of such gain scores seems to be ok:
"First, contrary to the traditional misconception, the reliability of gain scores is high in many practical situations, particularly when the pre- and posttest scores do not have equal variance and equal reliability."
so we have the following model:
library(tidyr) library(dplyr) data <- read.csv("http://pastebin.com/raw/G4D8dh1f") data_wide <- data %>% spread(time, undamaged_area) %>% separate(specimen_id, c("mat", "id", "tx")) %>% mutate(damage = before - after, unique_id = paste(mat, id, sep = "_")) %>% select(-mat, -tx, -id) # model for full factorial with replications mm2 <- lm(damage ~ material * treatment , data = data_wide)
The variance problem still remains. Confirmed by Levene's test:
library(car) leveneTest(damage ~ material * treatment, data_wide) # Levene's Test for Homogeneity of Variance (center = median) # Df F value Pr(>F) # group 9 4.8619 2.646e-05 *** # 90 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Following the link suggested by donlelek I found different approaches for anovas with heteroskedastic data. I tried to stabilize the variance by using log-transformation. Then Levene's test says that heterogeneity of the variance diappears:
data_wide <- within(data_wide, log_damage <- log(damage+1)) leveneTest(log_damage ~ material * treatment, data_wide) # Levene's Test for Homogeneity of Variance (center = median) # Df F value Pr(>F) # group 9 0.6916 0.7147 # 90 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The diagnostic plots seems not to be as weird as the previous ones (see donlelek's answer):
mm3 <- lm(log_damage ~ material * treatment, data_wide) plot(fitted(mm3), residuals(mm3, type = "pearson")) qqnorm(residuals(mm3, type = "pearson"))
The anova table gives the following output:
anova(mm3) # Analysis of Variance Table # # Response: log_damage # Df Sum Sq Mean Sq F value Pr(>F) # material 1 0.436 0.4362 0.7462 0.39 # treatment 4 83.652 20.9129 35.7786 < 2.2e-16 *** # material:treatment 4 20.213 5.0532 8.6452 5.966e-06 *** # Residuals 90 52.606 0.5845 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Just for double checking:
Anova function from
car package offers an option for heteroscedasticity correction. Interestingly, this function generates a roughly similar table for the "non-transformed"
Anova(mm2, white.adjust=TRUE) # Analysis of Deviance Table (Type II tests) # # Response: damage # Df F Pr(>F) # material 1 1.4251 0.2357 # treatment 4 28.2422 3.329e-15 *** # material:treatment 4 9.5739 1.701e-06 *** # Residuals 90 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
This double check gives me more confidence in the results. But do you think that the log-transformation is a reasonable approach? Can I trust the model now?