제가 맞춘 것은 다음과 같습니다:
Y ~ A(합계) + B(합계) + C(합계) + D(합계) + 주제(&랜덤)
JMP SLS(REML)와 Python statsmodels MixedLM(REML) 모두 동일한 데이터(~14,000개 관측치, ~200개 A 수준, 희소 A×C 교차표, ~120명의 무작위 피험자)를 사용하여 분석했습니다.
제가 관찰한 바는 다음과 같습니다.
- 요인 A에 대한 모든 약 200개의 LS 평균은 약 0.22만큼 일정하게 이동합니다(JMP 값이 더 높음).
- 상대적 효과와 순위는 동일합니다(차이 표준편차 < 0.001).
- JMP에서 다양한 요인에 대한 LS 평균은 다음과 같습니다. 서로 다른 큰 의미 (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, 윈도우.
안녕하세요 @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 )}}}
)
);
파이썬의 경우
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님 : 죄송하지만, 특히 파이썬 출력에서는 절편 등을 확인하기가 어렵습니다 (저는 파이썬을 사용하지 않습니다). JMP 출력 화면(절편, 계수, 최소제곱평균 등 포함)과 해당 파이썬 출력 화면을 나란히 배치하여 쉽게 비교할 수 있도록 보여주시겠습니까?
안녕하세요 @MRB3855 답변 감사합니다. 파이썬과 JMP 결과의 차이점을 시각화해 보았는데, 원인을 파악하는 데 도움이 될까요?
imageimageimage 영상
안녕하세요 @xmq님 : 처음에는 고정 효과에 대한 절편과 계수 간의 차이에 대한 모수화 문제일 수 있다고 생각했습니다( https://www.jmp.com/support/help/en/15.2/#page/jmp/nominal-factors.shtml# ).
하지만 LSMeans는 이와 무관해야 하므로 확실하지 않고, 저는 파이썬 사용자가 아니므로 다른 분이 도움을 주시면 좋겠습니다.
안녕하세요 @MRB3855 친절한 도움에 감사드립니다!
