我試穿了：
Y ~ A(Sum) + B(Sum) + C(Sum) + D(Sum) + Subject(&Random)
在 JMP SLS (REML) 和 Python statsmodels MixedLM (REML) 中，對同一資料集（~14k 個觀測值，~200 個 A 水平，稀疏 A×C 交叉表，~120 個隨機受試者）進行了分析。
我的觀察：
- 因子 A 的所有約 200 個最小平方法平均值均以約 0.22 的常數進行平移（JMP 較高）。
- 相對效應和排名完全相同（差異標準差 < 0.001）
- JMP中不同因子的最小平方法均值有 不同的大手段 （A → 78.5，B/D → 73.0），這在求和編碼中不應該發生。
- 同一個 X 座標下不同的 Y 座標基本上不會產生偏移（0.006）。
我需要幫助的地方：
JMP 的 SLS 模型（REML）為何會產生與標準 GLS β̂ = (X'V⁻¹X)⁻¹X'V⁻¹y 不同的截距？ SLS 模型在計算固定效應或最小平方法平均值時是否使用了與混合模型不同的演算法？
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 用戶，所以也許其他人可以幫忙。
