# Is it appropriate to use a multilevel model with little data?

Lets say I have v1 and v2 and I have 400 data points for each variable. I expect them to be negatively correlated, so a simple approach here would be to calculate the correlation

However, the 400 pairs of data points came from 40 different people, where each person only offered 10 pairs of datapoints

One way I could approach this would be to predict v1 from v2 and additionally treat person as a random effects variable, so something like v1 ~ v2 + (1|person)

I'm concerned about only have 10 datapoints per person. On the one hand a straight correlation doesnt consider the fact that the data is likely clustered due to multiple rows coming from the same individual. On the other hand, I only have 10 datapoints for an individual which may skew the correlation value for any given individual...

Is that a sufficient amount of data per level to run a multilevel model like this? Or does the N per level not matter, but rather the total N (400) being the only important factor?

Yes, it is reasonable to fit a multilevel model with this amount of data. Further, the single-level correlation you describe can be considered a special case of a multilevel model (one that assumes exactly zero person to person variance in the v1 intercept and the slope relating v1 to v2).

Before you fit this model, it may be worth asking: Do you expect the relationship between v1 and v2 to be driven by between-person variance, or by within-person variance? That is, does your theory suggest that people high in v1 will be low on v2, or that within a given person you will see instances such that where v1 is high, v2 will be low?

Perhaps your theory doesn't distinguish between these possibilities, but in principle there is no reason to suspect that the within-person relationship will be the same as the between-person relationship (i.e., see Simpson's Paradox and ecological fallacy).

To estimate both effects within the same model, you would compute means of v2 by person (v2pm), as well as instance-to-instance deflections from these means, or person-centered scores (v2pc).

data <- within(data, {v2pm = ave(v2, person, FUN=function(x) mean(x, na.rm=T))})

data\$v2pc <- data\$v2 - data\$v2pm You'd then fit a model such as: lmer(v1 ~ 1 + v2pm + v2pc + (1 + v2pc|person), data) If you primarily care about the effect of person-to-person differences in mean v2, you have 40 observations of these means. If you instead care about the effect of instance-to-instance differences in v2, you have 40 observations of such an effect (each based on 10 observations within a person). If you fit the model without decomposing the v2 person means from the person centered scores, the model will come to some weighted average of the between and within person effects. For more information on these kinds of models, see Bolger & Laurenceau 2013 and Gelman & Hill 2007. In a classical (frequentist/hypothesis testing) setting, this is IMHO a difficult tradeoff where neither of the solutions is better in all cases. Estimates with a random effect per person will be noisy, but calculating overall correlation looses information and may even be grossly misleading (e.g. Simpson's paradox). Which of the two is less wrong depends on your data and the decisions you want to use the data for. You basically want to take weighted average of the two, but I am not aware of any principled way to calculate the weights and determine p-values etc. (it definitely is possible, but I am not sure there are packages implementing it) In a fully Bayesian setting, hierarchical multilevel models can work well with little data per group. The "hierarchical" part implies something along the lines of:$X_i \sim N(\mu_i,\sigma)\mu_i = \beta Y_i + \alpha_{g(i)} \alpha \sim N(\nu, \tau)\nu \sim N(0,1)\tau \sim HalfNormal(0, 1)\beta \sim N(0,1)$Where$g(i)$indicates the group (person) the i-th datapoint belongs to. This is actually a principled approach to find a good compromise between overall and per-person correlation. This means that the values for$\alpha$are tied together via$\nu$(overall trend for all persons) and$\tau\$ (between person variance) and thus share information but are not forced to be equal. The scale parameters for priors (the ones) should be adjusted to reflect the scale of your data.

There is a worked out example going from classical linear regression to a hierarchical model in the Radon case study for Stan+Python. Also you may check out the documentation for R package rstanarm (especially the relevant vignette) that lets you fit these kind of models easily.

• However, the bit you did not address was few groups. That is in principle a potential problem. Of course, if you are willing to make assumption on the between person variability like a HN(0,1) prior, then this is less of an issue. – Björn Feb 14 '18 at 13:44
• I am not an expert, but the Stan case studies work with even fewer groups (e.g. the famous 8-schools example) and seem to be doing well. I guess few groups are dangerous when doing maximum a posterior or other point estimates, but when you perform full Bayes, the worst think that can happen is that you are (quite correctly) left with a large uncertainty in your posterior (and thus your inferences). – Martin Modrák Feb 14 '18 at 13:51
• thanks, but how about outside of a bayesian context? Im not familiar with bayes so that was never part of the analysis plan – Simon Feb 15 '18 at 0:31
• I edited the answer with my thoughts on the frequentist context. I think Bayes is the simplest way to take the structure of the data into account. Consider checking rstanarm out - it is IMHO easy to use (formula syntax as in R) and well documented. And don't be afraid of the word "Bayes" it just means that you use probabilities (not p-values,conf.int, Bayes factors, ...) to quantify uncertainty - as in "if our model is correct, the probability that the average effect is larger than 0.1 is 36%". – Martin Modrák Feb 15 '18 at 8:53