Executive summary
It is indeed often said that if all possible factor levels are included in a mixed model, then this factor should be treated as a fixed effect. This is not necessarily true FOR TWO DISTINCT REASONS:
(1) If the number of levels is large, then it can make sense to treat the [crossed] factor as random.
I agree with both @Tim and @RobertLong here: if a factor has a large number of levels that are all included in the model (such as e.g. all countries in the world; or all schools in a country; or maybe the entire population of subjects is surveyed, etc.), then there is nothing wrong with treating it as random --- this could be more parsimonious, could provide some shrinkage, etc.
lmer(size ~ age + subjectID) # fixed effect
lmer(size ~ age + (1|subjectID)) # random effect
(2) If the factor is nested within another random effect, then it has to be treated as random, independent of its number of levels.
There was a huge confusion in this thread (see comments) because other answers are about case #1 above, but the example you gave is an example of a different situation, namely this case #2. Here there are only two levels (i.e. not at all "a large number"!) and they do exhaust all possibilities, but they are nested inside another random effect, yielding a nested random effect.
lmer(size ~ age + (1|subject) + (1|subject:side) # side HAS to be random
Detailed discussion of your example
Sides and subjects in your imaginary experiment are related like classes and schools in the standard hierarchical model example. Perhaps each school (#1, #2, #3, etc.) has class A and class B, and these two classes are supposed to be roughly the same. You will not model classes A and B as a fixed effect with two levels; this would be a mistake. But you will not model classes A and B as a "separate" (i.e. crossed) random effect with two levels either; this would be a mistake too. Instead, you will model classes as a nested random effect inside schools.
See here: Crossed vs nested random effects: how do they differ and how are they specified correctly in lme4?
In your imaginary foot-size study, subject and side are random effects and side is nested inside subject. This essentially means that a combined variable is formed, e.g. John-Left, John-Right, Mary-Left, Mary-Right, etc., and there are two crossed random effects: subjects and subjects-sides. So for subject $i=1\ldots n$ and for side $j=1,2$ we would have:
$$\text{Size}_{ijk} = \mu+\alpha\cdot\text{Height}_{ijk}+\beta\cdot\text{Weight}_{ijk}+\gamma\cdot\text{Age}_{ijk}+\epsilon_i + \color{red}{\epsilon_{ij}} + \epsilon_{ijk}$$
$$\epsilon_i\sim\mathcal N(0,\sigma^2_\mathrm{subjects}),\quad\quad\text{Random intercept for each subject}$$
$$\color{red}{\epsilon_{ij}}\sim\mathcal N(0,\sigma^2_\text{subject-side}),\quad\quad\text{Random int. for side nested in subject}$$
$$\epsilon_{ijk}\sim\mathcal N(0,\sigma^2_\text{noise}),\quad\quad\text{Error term}$$
As you wrote yourself, "there is no reason to believe that right feet will on average be larger than left feet". So there should be no "global" effect (neither fixed nor random crossed) of right or left foot at all; instead, each subject can be thought of having "one" foot and "another" foot, and this variability we should include into the model. These "one" and "another" feet are nested within subjects, hence nested random effects.
More details in response to the comments. [Sep 26]
My model above includes Side as a nested random effect within Subjects. Here is an alternative model, suggested by @Robert, where Side is a fixed effect:
$$\text{Size}_{ijk} = \mu+\alpha\cdot\text{Height}_{ijk}+\beta\cdot\text{Weight}_{ijk}+\gamma\cdot\text{Age}_{ijk} + \color{red}{\delta\cdot\text{Side}_j}+\epsilon_i + \epsilon_{ijk}$$
I challenge @RobertLong or @gung to explain how this model can take care of the dependencies existing for consecutive measurements of the same Side of the same Subject, i.e. of the dependencies for data points with the same $ij$ combination.
It cannot.
The same is true for @gung's hypothetical model with Side as a crossed random effect:
$$\text{Size}_{ijk} = \mu+\alpha\cdot\text{Height}_{ijk}+\beta\cdot\text{Weight}_{ijk}+\gamma\cdot\text{Age}_{ijk} +\epsilon_i + \color{red}{\epsilon_j} + \epsilon_{ijk}$$
It fails to account for dependencies as well.
Demonstration via a simulation [Oct 2]
Here is a direct demonstration in R.
I generate a toy dataset with five subjects measured on both feet for five consecutive years. The effect of age is linear. Each subject has a random intercept. And each subject has one of the feet (either the left or the right) larger than another one.
set.seed(17)
demo = data.frame(expand.grid(age = 1:5,
side=c("Left", "Right"),
subject=c("Subject A", "Subject B", "Subject C", "Subject D", "Subject E")))
demo$size = 10 + demo$age + rnorm(nrow(demo))/3
for (s in unique(demo$subject)){
# adding a random intercept for each subject
demo[demo$subject==s,]$size = demo[demo$subject==s,]$size + rnorm(1)*10
# making the two feet of each subject different
for (l in unique(demo$side)){
demo[demo$subject==s & demo$side==l,]$size = demo[demo$subject==s & demo$side==l,]$size + rnorm(1)*7
}
}
plot(1:50, demo$size)
Apologies for my awful R skills. Here is how the data look like (each consecutive five dots is one feet of one person measured over the years; each consecutive ten dots are two feet of the same person):
Now we can fit a bunch of models:
require(lme4)
summary(lmer(size ~ age + side + (1|subject), demo))
summary(lmer(size ~ age + (1|side) + (1|subject), demo))
summary(lmer(size ~ age + (1|subject/side), demo))
All models include a fixed effect of age and a random effect of subject, but treat side differently.
Model 1: fixed effect of side. This is @Robert's model. Result: age comes out not significant ($t=1.8$), residual variance is huge (29.81).
Model 2: crossed random effect of side. This is @gung's "hypothetical" model from OP. Result: age comes out not significant ($t=1.4$), residual variance is huge (29.81).
Model 3: nested random effect of side. This is my model. Result: age is very significant ($t=37$, yes, thirty-seven), residual variance is tiny (0.07).
This clearly shows that side should be treated as a nested random effect.
Finally, in the comments @Robert suggested to include the global effect of side as a control variable. We can do it, while keeping the nested random effect:
summary(lmer(size ~ age + side + (1|subject/side), demo))
summary(lmer(size ~ age + (1|side) + (1|subject/side), demo))
These two models do not differe much from #3. Model 4 yields a tiny and insignificant fixed effect of side ($t=0.5$). Model 5 yields an estimate of side variance equal to exactly zero.
