How should I define my linear mixed model? Should I include frequency as fixed or random effect? In our study we measured a continuous outcome in multiple subjects, for both right and left ears, and at multiple frequencies (see image). We would like to examine the effect of sex and ear laterality on the continuous outcome. We decided to perform a linear mixed model as we would like to control for the correlation between ears (so we will include subject as random effect).
However, we know that the continuous outcome depends on the frequency (for example: at 1 kHz we expect a lower value than at 2 kHz. This association is not linear). Therefore, we have to take that into account.
I am not sure whether to include frequency as fixed effect or random effect? I am not interested in how frequency influences the outcome but I want to control for it? Could someone help me define the model in R?
I came up with the following models:
lme(outcome~ ear*gender*frequency, random= ~1|id, data=DF)

lme(outcome~ ear*gender, random= ~1|id, ~1|frequency, data=DF)


 A: Since you have 3500 subjects it is obvious that you need random intercepts for subject.
The larger question is whether you fit random intercepts for frequency.
It appears that subject and frequency are crossed - in that at least subjects were exposeed to different frequencies while for any particular frequency, multiple subjects were exposed to it. Note that your 2nd model it looks like you are attempting to fit random intercepts for both, however lme does not support crossed random effects and your model will be interpreted as specifying random intercepts for subject and a correlation structure based on frequency. To fit crossed random effects you would need to use lme4 or another package that supports crossed factors.
Since you manipulate the frequency that subjects are exposed to, this tends to point in the direction of modelling it as a fixed effect, not random, but since you know the association with your outcome is nonlinear you need to account for this. One way to do this, if you are unsure of the functional form of the association is to use a generalized additive mixed model. However, I would expect that you would have some idea of the functional form involved - so you might try to include relevant nonlinear terms or transforming the frequency variable. This of course assumes that you code frequency as continuous, not as a factor. If it is a factor then you can simply leave it as a fixed effect - but in that case you probably don't want to specify it in a three way interaction with ear and gender because that will result in a LOT of estimates and make interpretation difficult.
So my suggestion as a good place to start is with either:
lmer(outcome ~ ear*gender + frequency + (1|id), data=DF)

where frequency is either a factor, or if continuous then it has been transformed or includes nonlinear terms. If you believe that ear and gender and their interaction estimates change at different frequencies then you are going to need to include frequency as an interaction:
lmer(outcome ~ ear*gender*frequency + (1|id), data=DF)

As mentioned above, you probably want to avoid having frequency as a factor in this case. If something like a square root or log transform takes care of the nonlinearity, that should make the model interpretable. If you have to include higher order terms such as quadratic then you will need to decide whether to include those in the interaction as well.
Going back to the idea of treating frequency as random, then you would fit the model:
lmer(outcome ~ ear*gender + (1|id) + (1|frequencyID), data=DF)

however this will treat the "effect" of frequency as a normally distributed variate which may not sufficiently capture the nonlinearity.
I would try both approaches and if the answers to your research question(s) are broadly the same then all is good, and if they are different then you will have a very interesting problem to solve!
A: Your second model, lme(outcome~ ear*gender, random= ~1|id, ~1|frequency, data=DF), looks like the better model. If you are not interested in the effect of frequency, you still should control for frequency (perhaps as a factor since it is unlikely your outcome is linear in frequency).
Therefore, I would consider a model like the following:
lme(outcome~ ear*gender+frequency*gender, random= ~1|id, ~1|frequency, data=DF).
That way you will correct for frequency and frequency$\times$gender effects and your random effects will be centered around 0.
