# Cells means model in R for unbalanced Latin square design

I used a 4 x 4 Latin square experimental design recently and have been advised that a cell means - rather than sum of squares - approach is more appropriate because I have an unequal number of subjects in each cell.

I have found some information online which say I can fit a cell means model using the following:

mod <- lm(dv ~ 0+factor1:factor2, data = data)

Then is I use summary(mod), the coefficient estimates correspond to cell means of the Latin square.

There a few issues with this. First, I do not understand what the significance value of the t-statistic relates to and what I can infer from the significance value. Second, since this is a Latin square design, I only really want to test the main effects of the two blocking variables and treatment variable, as the test of the interaction will be biased.

So, my question is, how can I estimate a cell means model which includes the main effects of both my blocking and treatment variables and how should I then interpret the coefficients?

Any help or advice would be much appreciated. I can provide more details about the data and design.

• I think it would be helpful to describe in what way your design is unbalanced. By cell, do you mean a unique combination of rows and columns? If so, then, do you have multiple treatments in each cell? Or a single treatment in each cell with multiple observations? – Sal Mangiafico Oct 5 '17 at 21:42
• Yes, by cell I mean the unique combination of rows and columns. There is a single treatment per cell with multiple observations. – Con Des Oct 6 '17 at 0:45
• In your case, I don't see any problem with using the standard general linear model, model = lm(dv ~ treatment + row + column) . But you don't want to use the anova function. Use the Anova function in the car package. library(car); Anova(model) . – Sal Mangiafico Oct 6 '17 at 1:42
• It would be helpful if you could give the source or rationale for needing the cell means model. – Sal Mangiafico Oct 6 '17 at 1:44
• Hi @SalMangiafico - There is a chapter in Kirk (2014) on Latin square designs where he describes that this approach is appropriate when there are missing observations or unequal observations. I also confirmed this with him over email. Perhaps an additional complication is that in my design, treatment and column factors are within-subjects. Repeated-meaures Anova in R, I believe, only gives the type I SS and ezANOVA does not work for unbalanced Latin square designs. – Con Des Oct 6 '17 at 8:36

My recommendation would be to use a mixed effects model to account for the repeated measures nature of the design. I think the following may be a simple and viable solution. You may choose a more complicated model, or a statistician may find a different model more satisfying or appropriate.

These toy data are arranged in Latin square design (3 x 3), but with each subject subjected to 3 observations each. Each subject is given the same row (nominally random) for each observation, and then 3 random columns. This yields an unequal number of observation for each cell, and for each treatment.

Sources.

### Adapted from: http://rcompanion.org/handbook/I_07.html
###               http://rcompanion.org/handbook/C_04.html


Load the packages we'll need.

if(!require(lme4)){install.packages("lme4")}
if(!require(lmerTest)){install.packages("lmerTest")}
if(!require(multcompView)){install.packages("multcompView")}
if(!require(lsmeans)){install.packages("lsmeans")}
if(!require(ggplot2)){install.packages("ggplot2")}
if(!require(FSA)){install.packages("FSA")}


Then the data.

Input = ("
Obs Subject Row  Column  Rep  Measurement
1   a       R1   C1      1    3
2   a       R1   C2      2    14
3   a       R1   C3      3    19
4   b       R2   C1      1    11
5   b       R2   C3      2    11
6   b       R2   C2      3    17
7   c       R3   C2      1    8
8   c       R3   C2      2    9
9   c       R3   C1      3    18
10  d       R1   C3      1    17
11  d       R1   C2      2    14
12  d       R1   C3      3    19
13  e       R2   C1      1    11
14  e       R2   C1      2    9
15  e       R2   C1      3    16
16  f       R3   C2      1    8
17  f       R3   C2      2    9
18  f       R3   C2      3    13
19  g       R1   C3      1    15
20  g       R1   C3      2    13
21  g       R1   C3      3    17
22  h       R2   C1      1    11
23  h       R2   C2      2    16
24  h       R2   C3      3    12
25  i       R3   C1      1    20
26  i       R3   C2      2    10
27  i       R3   C3      3    15
")



I'll use a table of values to determine the treatment based on the row and column values, and then re-sort the data frame.

Rows = c(rep("R1", 3), rep("R2", 3), rep("R3", 3))
Columns = rep(c("C1", "C2", "C3"),3)
Treatment = c("A", "B", "C", "B", "C", "A", "C", "A", "B")
Square = data.frame(Rows, Columns, Treatment)

Data = merge(Data, Square, by.x=c("Row", "Column"), by.y=c("Rows", "Columns"))

Data = Data[order(Data\$Obs),]


Summarize the data. Note the number of observations in each cell.

library(FSA)

Summarize(Measurement ~ Row + Column, data = Data, digits=2)


The model is a mixed effects model. It assumes there are no temporal correlation effects within a subject. That is, Rep is not included in the model. But it does allow for there to be an effect of each subject. Subject is treated as a random variable. The mixed-effects model and type-iii sums-of-squares should handle the unbalanced observations okay.

library(lme4)

library(lmerTest)

model = lmer(Measurement ~ Treatment + Row + Column + (1|Subject),
data=Data,
REML=TRUE)

anova(model)

rand(model)


Least-square means are compared between each pair of treatments and summarized in a compact letter display. Least-square means can be more appropriate for unbalanced data than arithmetic means.

library(multcompView)

library(lsmeans)

leastsquare = lsmeans(model,
pairwise ~ Treatment,

CLD = cld(leastsquare,
alpha=0.05,
Letters=letters,

CLD


And a plot of the l.s. means and confidence intervals for treatments.

library(ggplot2)

qplot(x    = Treatment,
y    = lsmean,
data = CLD) +

geom_errorbar(aes(ymin  = lower.CL,
ymax  = upper.CL,
width = 0.15))

• I believe this fairly accurately reproduces my data - thank you! I had initially planned on using mixed effects but wasn't sure whether I would need to include a random effect for the column factor, as well. But I think this will work well. Thanks again for your help. – Con Des Oct 7 '17 at 10:26
• You can view your row or column variables as random variables. The way I usually think about this is: Let's say you were measuring effect of student learning for some curriculum, and you are measuring students at four different schools. You might include school as a blocking variable. If you care about the specific school ---- That is if you care about the effect of Springfield and Shelbyville per se ---- then that variable should be a fixed effect. But if the schools are just any four random schools, then that variable should be treated as a random variable. – Sal Mangiafico Oct 7 '17 at 12:35
• That's a very helpful way of explaining it. Since my column blocking variable is an event narrative that the participant reads, and there are any number of potential event narratives in the world, I guess it makes sense to include this as a random effect. – Con Des Oct 7 '17 at 13:23