装着しました：
Y ~ A(合計) + B(合計) + C(合計) + D(合計) + Subject(&ランダム)
同じデータ（約14,000の観測値、約200のAレベル、疎なA×Cクロス集計、約120のランダムな被験者）に対して、JMP SLS（REML）とPython statsmodels MixedLM（REML）の両方で分析しました。
私が観察したこと：
- 因子Aに関する約200個のLS平均値はすべて、定数約0.22（JMP値が高い）だけシフトされている。
- 相対的な効果と順位は同一である（差分標準偏差 < 0.001）
- JMPにおけるさまざまな要因のLS平均値は 異なるグランド平均 (A → 78.5、B/D → 73.0) は、Sum コーディングでは発生しないはずです。
- 同じX座標上でY座標が異なると、実質的にオフセットは発生しない（0.006）。
私が助けを必要としている点：
JMPのSLSモデル（REML）が、標準的なGLSモデルβ̂ = (X'V⁻¹X)⁻¹X'V⁻¹yとは異なる切片を生成する既知の理由はありますか？SLSは、固定効果またはLS平均を計算するために、混合モデルモデルとは異なるアルゴリズムを使用していますか？
JMP Pro 17、Windows版。
こんにちは、 @xmq さん、ようこそ！それぞれの関連する出力結果を見せていただけますか？そうすれば、何が起こっているのかを診断するのに役立ちます。いくつか思い当たる節はありますが、先走って結論を出すのは避けたいです。
こんにちは、 @MRB3855 さん、ご返信ありがとうございます！
ここに、私のデータセットとJMP/Pythonスクリプトおよび出力結果を共有します。
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 )}}}
)
);
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
ご意見やご提案があれば、ぜひお聞かせください。
こんにちは、@xmq さん。申し訳ありませんが、特に Python の出力からは切片などを確認するのが難しいです (私は Python を使用しません)。JMP の出力 (切片、係数、最小二乗平均などを含む) と、比較しやすいようにレイアウトされた対応する Python の出力の画像を見せていただけますか？
こんにちは、 @MRB3855さん。ご返信ありがとうございます。PythonとJMPの結果の違いを視覚化してみました。この違いの原因を特定するのに役立つでしょうか？
imageimageimage画像
こんにちは@xmq : 最初は、固定効果の切片と係数の差に関するパラメーター化の問題かもしれないと思いました ( https://www.jmp.com/support/help/en/15.2/#page/jmp/nominal-factors.shtml# )。
しかし、LSMeansはこの影響を受けないはずなので、確信が持てません。それに私はPythonユーザーではないので、どなたか助けていただけるとありがたいです。
こんにちは、 @MRB3855さん。親切なご協力ありがとうございます！
