946201_es-ES

liliconcom

¿Por qué SLS (REML) da una intersección diferente a la de statsmodels de Python para el mismo modelo mixto?

Me quedaba bien:

Y ~ A(Suma) + B(Suma) + C(Suma) + D(Suma) + Sujeto(&Aleatorio)

en JMP SLS (REML) y Python statsmodels MixedLM (REML) sobre los mismos datos (~14k observaciones, ~200 niveles de A, tabulación cruzada A×C dispersa, ~120 sujetos aleatorios).

Lo que observo:

• Todas las ~200 medias de mínimos cuadrados para el factor A se desplazan en una constante de ~0,22 (JMP más alto).
• Los efectos relativos y las clasificaciones son idénticos (dist. < 0,001).
• Las medias de mínimos cuadrados para diferentes factores en JMP tienen diferentes grandes medios (A → 78.5, B/D → 73.0), lo cual no debería ocurrir con la codificación Sum.
• Un valor diferente de Y en el mismo X no produce prácticamente ningún desplazamiento (0,006).

En qué necesito ayuda:

¿Existe alguna razón conocida por la que la configuración SLS (REML) de JMP produzca una intersección diferente a la del GLS estándar β̂ = (X'V⁻¹X)⁻¹X'V⁻¹y? ¿Utiliza SLS un algoritmo diferente al de la configuración de Modelo Mixto para calcular los efectos fijos o las medias de mínimos cuadrados?

JMP Pro 17, Windows.

Modelado estadístico avanzadoRe: ¿Por qué SLS (REML) da una intersección diferente a la de statsmodels de Python para el mismo modelo mixto?

Hola @xmq , ¡bienvenido! ¿Podrías mostrar la información relevante de cada comando, por favor? Nos ayudará a diagnosticar el problema. Tengo algunas ideas, pero no quiero adelantarnos a los acontecimientos.

Re: ¿Por qué SLS (REML) da una intersección diferente a la de statsmodels de Python para el mismo modelo mixto?

Hola @MRB3855 , ¡gracias por tu respuesta!

Aquí comparto mi conjunto de datos y los scripts y resultados de JMP/Python para su verificación.

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

Para 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

Agradecería cualquier comentario o sugerencia que pudieran tener.

Re: ¿Por qué SLS (REML) da una intersección diferente a la de statsmodels de Python para el mismo modelo mixto?

Hola @xmq : Disculpa, pero me resulta difícil ver la intersección, etc., especialmente en la salida de Python (no uso Python). ¿Podrías mostrar una imagen de la salida de JMP (incluyendo la intersección, los coeficientes, las medias de mínimos cuadrados, etc.) y una imagen de la salida de Python correspondiente, organizadas de forma que sean fácilmente comparables?

Re: ¿Por qué SLS (REML) da una intersección diferente a la de statsmodels de Python para el mismo modelo mixto?

Hola @MRB3855 Gracias por tu respuesta, aquí he visualizado la diferencia entre los resultados de Python y JMP, ¿sería útil identificar la razón detrás de esto?
imageimageimage imagen

Re: ¿Por qué SLS (REML) da una intersección diferente a la de statsmodels de Python para el mismo modelo mixto?

Hola @xmq : al principio pensé que podría ser un problema de parametrización para las diferencias entre las intersecciones y los coeficientes para los efectos fijos ( https://www.jmp.com/support/help/en/15.2/#page/jmp/nominal-factors.shtml# ).

Pero no estoy tan seguro, ya que LSMeans no debería verse afectado por esto, y no soy usuario de Python, así que tal vez alguien más podría ayudar aquí.

Re: ¿Por qué SLS (REML) da una intersección diferente a la de statsmodels de Python para el mismo modelo mixto?

Hola @MRB3855 ¡Gracias por tu amable ayuda!