2
$\begingroup$

I've been looking into using mixed effect models for a project, and along the way wanted to experiment with how they'd work on synthetic data.

However, I'm confused by the results — if I generate what I think is null data (e.g. no true difference between two inputs), the p-value for their difference is skewed towards higher significance than I would expect.

I'm not sure if I'm making an incorrect assumption about how mixed effect models are supposed to work, or if I'm generating it incorrectly. Here's what I've been doing:

I have two sets of inputs, each with N distinct groups, and M samples per group. To draw outputs, for both sets of inputs, I draw N values to represent each group's baseline effect, and then M values per group to represent the random variances. I then add those up for both sets of inputs, run a mixed effect model, and generate a p-value. I repeated this 1000 times and made a histogram of the p-values.

Here's the code:

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

num_groups = 10
num_samples_per_group = 20
num_iter = 1000

group_var = 1.0
random_var = 1.0

# Build basic dataframe w/ two inputs (input_0 and input_1), each with num_groups different groups, each
# containing num_samples_per_group individual samples.
df = pd.DataFrame({
    'input': ['input_0'] * num_samples_per_group * num_groups + ['input_1'] * num_samples_per_group * num_groups,
    'group': (
        ['input_0_group_%s' % (i % num_samples_per_group) for i in range(num_samples_per_group * num_groups)] +
        ['input_1_group_%s' % (i % num_samples_per_group) for i in range(num_samples_per_group * num_groups)])})

pvals = []
for i in range(num_iter):
    # Draw a group effect for each group w/ variance group_var
    input_one_group_effects = np.random.normal(0, np.sqrt(group_var), num_groups)
    input_two_group_effects = np.random.normal(0, np.sqrt(group_var), num_groups)

    # Draw a random effect for each sample in each group w/ variance random_var
    input_one_random_effects = np.random.normal(0, np.sqrt(random_var), (num_samples_per_group, num_groups))
    input_two_random_effects = np.random.normal(0, np.sqrt(random_var), (num_samples_per_group, num_groups))
    
    # Add group and random effects.
    input_one_effects = (input_one_group_effects + input_one_random_effects).reshape((num_samples_per_group * num_groups,))
    input_two_effects = (input_two_group_effects + input_two_random_effects).reshape((num_samples_per_group * num_groups,))
    
    # Insert into dataframe and run stats model.
    df['output'] = np.concatenate((input_one_effects, input_two_effects))
    model = sm.MixedLM.from_formula(
            "output ~ input",
            data=df,
            groups="group",
    ).fit()
    pvals.append(model.pvalues[1])

plt.hist(pvals)

And here's the resulting histogram:

Histogram of p-values from the above snippet

I'm not sure if this is because I'm misunderstanding the assumptions of mixed effect models, or if this is expected when there's limited #'s of groups?

Thanks a bunch for your help.

$\endgroup$
0
$\begingroup$

The situation with a small number of groups is very challenging. If you switch it to 20 groups and 10 samples per group (instead of 10 groups and 20 samples per group), the performance will be much better.

Right now we use a normal reference distribution to get the p-values. We should add an option to use the t-distribution. There is a literature on degrees of freedom for this (Satterthwaite, Kenward-Roger, etc.), which we have not implemented yet.

$\endgroup$
2
  • $\begingroup$ Thanks, that's very helpful -- does this mean that in a real-world situation where we want to model a similar setup to assess significance, we shouldn't rely on p-values below a certain # of groups? Or is there an alternative you'd suggest? $\endgroup$ – Quat Cola Jan 22 at 18:36
  • $\begingroup$ The performance depends not only on the number of groups but also on the variance (group and individual). The Satterthwaite and Kenward-Roger approaches are the best way to handle this, but they are tricky to implement (some details in appendix A here: jstatsoft.org/article/view/v082i13/0). I will get to it eventually... With few groups you could consider fixed effects for groups. Another alternative is GEE, which also has problems with a small number of large groups (see the work of Mancl, Westgate, etc.). $\endgroup$ – Kerby Shedden Jan 23 at 21:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.