946201_it-IT

liliconcom

Perché SLS (REML) fornisce un'intercetta diversa rispetto a statsmodels di Python per lo stesso modello misto?

Mi sono adattato:

Y ~ A(Somma) + B(Somma) + C(Somma) + D(Somma) + Soggetto(&Casuale)

in entrambi i modelli JMP SLS (REML) e Python statsmodels MixedLM (REML) sugli stessi dati (~14k osservazioni, ~200 livelli di A, tabulazione incrociata A×C sparsa, ~120 soggetti casuali).

Ciò che osservo:

• Tutte le ~200 medie LS per il fattore A sono traslate di una costante ~0,22 (JMP più alto)
• Gli effetti relativi e le classifiche sono identici (diff std < 0,001)
• I valori medi dei minimi quadrati per diversi fattori in JMP hanno diversi grandi mezzi (A → 78,5, B/D → 73,0), cosa che non dovrebbe accadere con la codifica della somma
• Un diverso valore di Y sullo stesso valore di X non produce praticamente alcuno scostamento (0,006)

Ho bisogno di aiuto per:

Esiste una ragione nota per cui la personalità SLS di JMP (REML) produce un'intercetta diversa rispetto al GLS standard β̂ = (X'V⁻¹X)⁻¹X'V⁻¹y? SLS utilizza un algoritmo diverso rispetto alla personalità del modello misto per il calcolo degli effetti fissi o delle medie LS?

JMP Pro 17, Windows.

Modellazione statistica avanzataOggetto: Perché SLS (REML) fornisce un'intercetta diversa rispetto a statsmodels di Python per lo stesso modello misto?

Ciao @xmq e benvenuto! Potresti mostrarci l'output pertinente di ciascuno, per favore? Aiuterebbe gli altri a capire cosa potrebbe succedere; ho qualche idea, ma non vorrei anticipare i tempi.

Oggetto: Perché SLS (REML) fornisce un'intercetta diversa rispetto a statsmodels di Python per lo stesso modello misto?

Ciao @MRB3855 , grazie per la risposta!

Qui condivido il mio dataset e gli script JMP/Python e gli output per la verifica,

Per 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 )}}}
	)
);

Per 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

Sarei grato per qualsiasi commento o suggerimento.

Oggetto: Perché SLS (REML) fornisce un'intercetta diversa rispetto a statsmodels di Python per lo stesso modello misto?

Ciao @xmq : Mi scuso, ma mi è difficile visualizzare l'intercetta ecc. soprattutto dall'output di Python (non uso Python). Potresti mostrarmi un'immagine dell'output di JMP (inclusi intercetta, coefficienti, ls means, ecc.) e un'immagine del corrispondente output di Python disposti in modo che siano facilmente confrontabili?

Oggetto: Perché SLS (REML) fornisce un'intercetta diversa rispetto a statsmodels di Python per lo stesso modello misto?

Ciao @MRB3855 , grazie per la risposta. Qui ho visualizzato la differenza tra i risultati di Python e JMP. Sarebbe utile capire il motivo di questa differenza?
imageimageimage immagine

Oggetto: Perché SLS (REML) fornisce un'intercetta diversa rispetto a statsmodels di Python per lo stesso modello misto?

Ciao @xmq : inizialmente ho pensato che potesse trattarsi di un problema di parametrizzazione per le differenze tra le intercette e i coefficienti per gli effetti fissi ( https://www.jmp.com/support/help/en/15.2/#page/jmp/nominal-factors.shtml# ).

Però non ne sono così sicuro, dato che gli LSMean non dovrebbero essere influenzati da questo, e non sono un utente Python, quindi forse qualcun altro potrebbe aiutarmi.

Oggetto: Perché SLS (REML) fornisce un'intercetta diversa rispetto a statsmodels di Python per lo stesso modello misto?

Ciao @MRB3855, grazie per il tuo gentile aiuto!