Fixed effect vs random effect when all possibilities are included in a mixed effects model In a mixed effects model the recommendation is to use a fixed effect to estimate a parameter if all possible levels are included (e.g., both males and females).  It is further recommended to use a random effect to account for a variable if the levels included are just a random sample from a population (enrolled patients from the universe of possible patients) and you want to estimate the population mean and variance instead of the means of the individual factor levels.  
I am wondering if you are logically obliged to always use a fixed effect in this manner.  Consider a study on how foot / shoe size changes through development and is related to, say, height, weight and age.  ${\rm Side}$ clearly must be included in the model somehow to account for the fact that the measurements over the years are nested within a given foot and are not independent.  Moreover, right and left are all the possibilities that can exist.  In addition, it can be very true that for a given participant their right foot is larger (or smaller) than their left.  However, although foot size does differ somewhat between the feet for all people, there is no reason to believe that right feet will on average be larger than left feet.  If they are in your sample, this is presumably due to something about the genetics of the people in your sample, rather than something intrinsic to right-foot-ness.  Finally, ${\rm side}$ seems like a nuisance parameter, not something you really care about.  
Let me note that I made this example up.  It may not be any good; it is just to get the idea across.  For all I know, having a large right foot and a small left foot was necessary for survival in the paleolithic. 
In a case like this, would it make (more / less / any) sense to incorporate ${\rm side}$ in the model as a random effect?  What would be the pros and cons of using a fixed vs. random effect here?  
 A: To add to the other answers:
I don't think you are logically obliged to always use a fixed effect in the manner described in the OP. Even when the usual definitions/guidelines for when to treat a factor as random are not met, I might be inclined to still model it as random when there are a large number of levels, so that treating the factor as fixed would consume many degrees of freedom and result in a cumbersome and less parsimonious model. 
A: If you're talking about the situation where you know all possible levels of a factor of interest, and also have data to estimate the effects, then definitely you don't need to represent levels with random effects. 
The reason that you want to set random effect to a factor is because you wish to make inference on the effects of all levels of that factor, which are typically unknown. To make that kind of inference, you impose the assumption that the effects of all levels form a normal distribution in general. But given your problem setting, you can estimates the effects of all levels. Then there is certainly no need to set random effects and impose additional assumption. 
It's like the situation that you are able to get all the values of the population (thus you know the true mean), but you are trying to take a large sample from the population and use central limit theorem to approximate the sampling distribution, and then make inference on the true mean.
A: 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.
A: The general problem with "fixed" and "random" effects is that they are not defined in a consistent way. Andrew Gelman quotes several of them:

(1) Fixed effects are constant across individuals, and random effects
  vary. For example, in a growth study, a model with random intercepts
  $a_i$ and fixed slope $b$ corresponds to parallel lines for different
  individuals $i$, or the model $y_{it} = a_i + b_t$. Kreft and De Leeuw
  (1998) thus distinguish between fixed and random coefficients.
(2) Effects are fixed if they are interesting in themselves or random
  if there is interest in the underlying population. Searle, Casella,
  and McCulloch (1992, Section 1.4) explore this distinction in depth.
(3) “When a sample exhausts the population, the corresponding variable
  is fixed; when the sample is a small (i.e., negligible) part of the
  population the corresponding variable is random.” (Green and Tukey,
  1960)
(4) “If an effect is assumed to be a realized value of a random
  variable, it is called a random effect.” (LaMotte, 1983)
(5) Fixed effects are estimated using least squares (or, more
  generally, maximum likelihood) and random effects are estimated with
  shrinkage (“linear unbiased prediction” in the terminology of
  Robinson, 1991). This definition is standard in the multilevel
  modeling literature (see, for example, Snijders and Bosker, 1999,
  Section 4.2) and in econometrics.

and notices that they are not consistent. In his book Data Analysis Using Regression and Multilevel/Hierarchical Models he generally avoids using those terms and in their work he focuses on fixed or varying between groups intercepts and slopes because

Fixed effects can be viewed as special cases of random effects, in
  which the higher-level variance (in model (1.1), this would be
  $\sigma^2_\alpha$ ) is set to $0$ or $\infty$. Hence, in our framework,
  all regression parameters are “random,” and the term “multilevel” is
  all-encompassing.

This is especially true with Bayesian framework - commonly used for mixed models - where all the effects are random per se. If you are thinking Bayesian, you are not really concerned with "fixed" effects and point estimates and have no problem with treating all the effects as random.
The more I read on this topic, the more I am convinced that this is rather an ideological discussion on what we can (or should) estimate and what we only can predict (here I could refer also to your own answer). You use random effects if you have a random sample of possible outcomes, so you are not concerned about individual estimates and you care rather about the population effects, then individuals. So the answer of your question depends also on what do you think about if you want or can estimate the fixed effects given your data. If all the possible levels are included in your data you can estimate fixed effects - also, like in your example, the number of levels could be small and that would generally not be good for estimating random effects and there are some minimal requirements for this.
Best case scenario argument
Say you have unlimited amounts of data and unlimited computational power. In this case you could imagine estimating every effect as fixed, since fixed effects give you more flexibility (enable us to compare the individual effects). However, even in this case, most of us would be reluctant to use fixed effects for everything.
For example, imagine that you want to model exam results of schools in some region and you have data on all the 100 schools in the region. In this case you could threat schools as fixed - since you have data on all the levels - but in practice you probably would rather think of them as random. Why is that?


*

*One reason is that generally in this kind of cases you are not interested in effects of individual schools (and it is hard to compare all of them), but rather a general variability between schools.

*Another argument in here is model parsimony. Generally you are not interested in "every possible influence" model, so in your model you include few fixed effects that you want to test and control for the other possible sources of variability. This makes mixed effects models fit the general way of thinking about statistical modeling where you estimate something and control for other things. With complicated (multilevel or hierarchical) data you have many effects to include, so you threat some as "fixed" and some as "random" so to control for them.

*In this scenario you also wouldn't think of the schools as each having its own, unique, influence on the results, but rather as about schools having some influence in general. So this argument would be that we believe that is is not really possible to estimate the unique effects of individual schools and so we threat them as random sample of possible schools effects.
Mixed effects models are somewhere in between "everything fixed" and "everything random" scenarios. The data we encounter makes us to lower our expectations about estimate everything as fixed effects, so we decide what effects we want to compare and what effects we want to control, or have general feeling about their influence. It is not only about what the data is, but also how we think of the data while modeling it.
A: Following the above discussion, I thought that side could also be modelled on a random slope across subjects, that is with the following LME model:
lmer(size ~ age + (1+side|subject), demo) [model 4 or lme4]
(because as you said: there is a random variation of size across subject and in addition to this, a random variation of the side effect across subjects).
I checked this model in the simulation above. It gives the same result as model 3 (noted lme3 below) (where side was modelled as a random factor nested in subject): t for Age = 37, residual variance = 0.07. The dotplot helps understanding the similarity between the two models :

*

*For model 3: dotplot(ranef(lme3, condVar=T))$side



*

*For model 4: dotplot(ranef(lme4, condVar=T))


I found this very enlightening and thought I would share it. This is my first participation here, so I hope I am not missing the point. best, N
