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:

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.