1
$\begingroup$

I'm trying to analyze a data set but I'm having trouble finding the right analysis. I'm looking at performance (Work) of a material.

Unfortunately my experimental design turned out to be unbalanced. I had 16 individuals set up in one of two treatments (A or B) (between treatments). After treatment we collected material and switch the individuals to the other treatment and collected material again. Individuals were switch about 6 times, and I was not able to collect material each time.

After collection I measured performance (Work) at four different conditions (a, b, c or d)(within treatment).

I want to see if Treatment had an effect on the performance of the material in any of the conditions. So, I ran a linear mixed model on the log transformed Work data using R. The original Work data did not met ANOVA assumptions. My response variable is Log of Work, my factors are between treatment and within treatment, and I used individual and order as my random effects.

I set zero contrast to factors, including random effects. Is that OK? Should I use a different analysis for my type of data?

Thank you for your help.

Here is the data file: https://docs.google.com/spreadsheets/d/1NnYWbvpjoCLAQJNaGnFDEf_DHdDPtRgjHSZdSe6ChVI/edit?usp=sharing

Here is my code:

library(lme4) # for glmer
library(multcomp) # for glht
library(emmeans) # for emm

# read in a data file 
Mchnqs = read.csv("Example.csv")
colnames(Mchnqs)[1] <- gsub('^...','',colnames(Mchnqs)[1])
View(Mchnqs)
Mchnqs$Individual = factor(Mchnqs$Individual) # convert to nominal factor
Mchnqs$Order = factor(Mchnqs$Order) # convert to nominal factor
Mchnqs$Within.Treatment = factor(Mchnqs$Within.Treatment) # convert to nominal factor
Mchnqs$Betweenn.Treatment = factor(Mchnqs$Betweenn.Treatment) # convert to nominal factor
summary(Mchnqs)
with(Mchnqs, interaction.plot(Betweenn.Treatment, Within.Treatment, Log.Work)) # for convenience
boxplot(Log.Work ~ Betweenn.Treatment*Within.Treatment, data=Mchnqs, xlab="Between Treatment", ylab="Work (Log uJ)") # boxplots


# set sum-to-zero contrasts 
contrasts(Mchnqs$Betweenn.Treatment) <- "contr.sum"
contrasts(Mchnqs$Within.Treatment) <- "contr.sum"
contrasts(Mchnqs$Order) <- "contr.sum"
contrasts(Mchnqs$Individual) <- "contr.sum"

# LMM test on Log.Work
# Betweenn.Treatment:Within.Treatment are
# fixed effects. Individual and Order are random effects.
m = lmer(Log.Work ~ (Betweenn.Treatment * Within.Treatment) + (1|Order) + (1|Individual), data=Mchnqs)
Anova(m, type=3)
qqnorm(residuals(m)); qqline(residuals(m)) # plot residuals

# perform post hoc pairwise comparisons with holm's bonferroni correction
summary(glht(m, emm(pairwise ~ Betweenn.Treatment * Within.Treatment)), test=adjusted(type="holm"))
$\endgroup$

1 Answer 1

3
$\begingroup$

Linear mixed models work well with unbalanced data, unlike ANOVA. Note that the assumptions of lmms are different to ANOVA meaning that we do NOT assume normality or homoscedasticity of the response variable. Lmms make, f.e., assumptions on the errors. My suggestion is to fit the model again with an untransformed response variable, then check your assumptions graphically with a qqplot. See this thread Checking assumptions lmer/lme mixed models in R, for example.

I am not sure whether it makes sense to fit a random intercept for order, maybe a random slope is more suitable, as this is a within-design.

I also suggest that if you would like to do ANOVA with an unbalanced design, first use type III sums of squares, if an interaction is present stop, if not, test main effects again with type II.

$\endgroup$
2
  • $\begingroup$ Thank you for your answer. $\endgroup$ Commented Mar 7, 2022 at 18:24
  • $\begingroup$ Thank you for your answer. The qq plot of the residuals for the untrasformed after the lmer function shows a positive skewed distribution. That's partly why I used log-transformed data. I think I could drop the order random variable from the model because I think it represents the same information as the individual random variable. I did not see an interaction effect. When I ran a type 2 ANOVA it gave me the same result as with a type 3. $\endgroup$ Commented Mar 7, 2022 at 18:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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