Ich passte:
Y ~ A(Summe) + B(Summe) + C(Summe) + D(Summe) + Subjekt(&Zufällig)
Sowohl in JMP SLS (REML) als auch in Python statsmodels MixedLM (REML) auf den gleichen Daten (~14k Beobachtungen, ~200 Stufen von A, spärliche A×C-Kreuztabelle, ~120 zufällige Probanden).
Was ich beobachte:
- Alle ~200 LS-Mittelwerte für Faktor A sind um eine Konstante von ~0,22 verschoben (JMP höher).
- Die relativen Effekte und Rangfolgen sind identisch (Differenz Standardabweichung < 0,001).
- LS-Mittelwerte für verschiedene Faktoren in JMP haben verschiedene große Mittel (A → 78,5, B/D → 73,0), was bei der Summenkodierung nicht vorkommen sollte.
- Ein anderes Y bei gleichem X ergibt praktisch keinen Offset (0,006).
Wobei ich Hilfe benötige:
Gibt es einen bekannten Grund, warum die SLS-Persönlichkeit von JMP (REML) einen anderen Achsenabschnitt liefert als die Standard-GLS β̂ = (X'V⁻¹X)⁻¹X'V⁻¹y? Verwendet SLS einen anderen Algorithmus als die Mixed-Model-Persönlichkeit zur Berechnung der festen Effekte oder der LS-Mittelwerte?
JMP Pro 17, Windows.
Hallo @xmq und herzlich willkommen! Könntest du bitte die relevanten Ausgaben der einzelnen Befehle zeigen? Das würde uns hier helfen, die Ursache des Problems zu finden. Ich habe zwar schon ein paar Vermutungen, aber ich möchte nicht voreilig Schlüsse ziehen.
Hallo @MRB3855 , vielen Dank für deine Antwort!
Hier stelle ich meinen Datensatz sowie meine JMP/Python-Skripte und -Ausgaben zur Überprüfung zur Verfügung.
Für JMP
Fit Model(
Y( :Y1 ),
Effects( :A, :B, :C, :D ),
Random Effects( :Subject ),
NoBounds( 1 ),
Personality( "Standard Least Squares" ),
Method( "REML" ),
Emphasis( "Effect Leverage" ),
Run(
:Y1 << {Summary of Fit( 1 ), Analysis of Variance( 0 ),
Parameter Estimates( 1 ), Scaled Estimates( 0 ),
Plot Actual by Predicted( 0 ), Plot Regression( 0 ),
Plot Residual by Predicted( 0 ), Plot Studentized Residuals( 0 ),
Plot Effect Leverage( 0 ), Plot Residual by Normal Quantiles( 0 ),
{:A << {LSMeans Student's t( 0.05 )}}}
)
);
Für Python
def fit_model(df, y_col):
"""Fit Y ~ A(Sum) + B(Sum) + C(Sum) + D(Sum) + (1|Subject) via REML."""
# Note: column "C" conflicts with patsy's C() function,
# so we temporarily rename it for formula construction.
df_tmp = df.rename(columns={'C': '_C_'})
formula = '1 + C(A, Sum) + C(B, Sum) + C(_C_, Sum) + C(D, Sum)'
X_mat = dmatrix(formula, df_tmp)
design_info = X_mat.design_info
X = pd.DataFrame(np.asarray(X_mat), columns=design_info.column_names)
y = df[y_col].values
groups = df['Subject']
model = MixedLM(y, X, groups=groups)
result = model.fit(reml=True, method='powell')
return result, X, design_info, df_tmp
def compute_lsmeans_A(result, X, df_tmp, design_info):
"""Compute LS Means for factor A (Sum coding → intercept + effect)."""
fe = result.fe_params.values
n_fe = len(fe)
col_names = list(X.columns)
target_prefix = 'C(A, Sum)'
target_col_idx = [i for i, c in enumerate(col_names) if target_prefix in c]
base_contrast = np.zeros(n_fe)
base_contrast[0] = 1.0 # intercept
levels = sorted(df_tmp['A'].unique())
lsmeans = {}
for level in levels:
sr = df_tmp.iloc[0:1].copy()
sr['A'] = level
X_single = np.asarray(build_design_matrices([design_info], sr)[0])
c = base_contrast.copy()
for idx in target_col_idx:
c[idx] = X_single[0, idx]
lsmeans[level] = float(c @ fe)
return lsmeans
Ich freue mich über Ihre Kommentare/Vorschläge.
Hallo @xmq , Entschuldigung, aber es fällt mir schwer, den Achsenabschnitt usw., insbesondere in der Python-Ausgabe, zu erkennen (ich verwende kein Python). Könnten Sie bitte ein Bild der JMP-Ausgabe (einschließlich Achsenabschnitt, Koeffizienten, Mittelwert usw.) und ein Bild der entsprechenden Python-Ausgabe zeigen, sodass sie gut vergleichbar sind?
Hallo @MRB3855, vielen Dank für Ihre Antwort. Hier habe ich den Unterschied zwischen den Ergebnissen von Python und JMP visualisiert. Wäre es hilfreich, die Ursache dafür zu ermitteln?
imageimageimage Bild
Hallo @xmq : Zuerst dachte ich, es könnte ein Problem mit der Parametrisierung sein, da es Unterschiede zwischen den Achsenabschnitten und Koeffizienten für die festen Effekte gibt ( https://www.jmp.com/support/help/en/15.2/#page/jmp/nominal-factors.shtml# ).
Aber ich bin mir da nicht so sicher, da die LSMeans davon unberührt bleiben sollten, und ich bin kein Python-Anwender, daher kann vielleicht jemand anderes hier weiterhelfen.
Hallo @MRB3855, vielen Dank für Ihre freundliche Hilfe!
