J'ai ajusté :
Y ~ A(Somme) + B(Somme) + C(Somme) + D(Somme) + Sujet(&Aléatoire)
dans JMP SLS (REML) et Python statsmodels MixedLM (REML) sur les mêmes données (~14k obs, ~200 niveaux de A, tableau croisé A×C clairsemé, ~120 sujets aléatoires).
Ce que j'observe :
- Les ~200 moyennes LS pour le facteur A sont toutes décalées d'une constante d'environ 0,22 (JMP plus élevé).
- Les effets relatifs et les classements sont identiques (différence standard < 0,001).
- Les moyennes LS pour différents facteurs dans JMP ont différents grands moyens (A → 78,5, B/D → 73,0), ce qui ne devrait pas se produire avec le codage par somme.
- Une valeur Y différente sur la même valeur X ne donne pratiquement aucun décalage (0,006).
J'ai besoin d'aide pour :
Existe-t-il une raison connue pour laquelle le modèle SLS (REML) de JMP produit une ordonnée à l'origine différente de celle du modèle GLS standard β̂ = (X'V⁻¹X)⁻¹X'V⁻¹y ? Le modèle SLS utilise-t-il un algorithme différent du modèle mixte pour le calcul des effets fixes ou des moyennes des moindres carrés ?
JMP Pro 17, Windows.
Salut @xmq , et bienvenue ! Pourriez-vous afficher les résultats pertinents de chaque commande ? Cela aidera l'équipe à identifier le problème. J'ai quelques pistes, mais je préfère ne pas tirer de conclusions hâtives.
Salut @MRB3855 , merci pour ta réponse !
Je partage ici mon jeu de données ainsi que mes scripts JMP/Python et leurs résultats pour vérification.
Pour 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 )}}}
)
);
Pour 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
Vos commentaires et suggestions seraient les bienvenus.
Bonjour @xmq : Veuillez m'excuser, mais il m'est difficile de visualiser l'ordonnée à l'origine, etc., notamment à partir de la sortie Python (je n'utilise pas Python). Pourriez-vous me fournir une capture d'écran de la sortie JMP (avec l'ordonnée à l'origine, les coefficients, les moyennes des moindres carrés, etc.) et une capture d'écran de la sortie Python correspondante, présentées de manière à ce qu'elles soient facilement comparables ?
Bonjour @MRB3855 , merci pour votre réponse. J'ai visualisé ici la différence entre les résultats de Python et de JMP. Serait-il utile d'en identifier la raison ?
imageimageimage image
Bonjour @xmq : au début, je pensais qu'il pouvait s'agir d'un problème de paramétrisation pour les différences entre les ordonnées à l'origine et les coefficients des effets fixes ( https://www.jmp.com/support/help/en/15.2/#page/jmp/nominal-factors.shtml# ).
Mais je n'en suis pas si sûr, car les LSMeans ne devraient pas être affectés par cela, et je n'utilise pas Python, donc peut-être que quelqu'un d'autre pourrait m'aider.
Bonjour @MRB3855, merci beaucoup pour votre aide précieuse !
