import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.regression.mixed_linear_model import MixedLM


# Standardize movement_time for easier interpretation 
# df contains columns: subject_ID, hemisphere, cortical_label,decoder_accuracy, movement_time

# (coefficient will be per 1 SD increase in movement_time)
df['movement_time_std'] = (df['movement_time'] - df['movement_time'].mean()) / df['movement_time'].std()


print("\n" + "="*80)
print("MODEL 2: With Label × Side Interactions")
print("Formula: accuracy ~ movement_time_std * (Label + side) + (1 | subject_ID)")
print("="*80)

model = smf.mixedlm("accuracy ~ movement_time_std * (Label + side)", 
                     data=df, 
                     groups=df["subject_ID"],
                     re_formula="~ 1")
result = model.fit(method='powell')

print(result.summary())

print("\nMODEL 2 KEY RESULTS:")
print(f"Fixed Effect (movement_time_std): {result.fe_params['movement_time_std']:.4f}")
print(f"95% CI: [{result.conf_int().loc['movement_time_std', 0]:.4f}, {result.conf_int().loc['movement_time_std', 1]:.4f}]")
print(f"Log-likelihood: {result.llf:.2f}")