General Linear Model¶

This example shows how to simulate and recover parameters for a simple linear regression model. Each subject has their own regression coefficients and we fit all subjects simultaneously using EMModel.

In [1]:

import numpy as np
from pyem import EMModel
from pyem.utils import plotting
from pyem.models.glm import glm_model # GLM

print(f"Model ID   : {glm_model.id}")
print(f"Description: {glm_model.desc}")
print(f"Spec       : {glm_model.spec}")

import numpy as np from pyem import EMModel from pyem.utils import plotting from pyem.models.glm import glm_model # GLM print(f"Model ID : {glm_model.id}") print(f"Description: {glm_model.desc}") print(f"Spec : {glm_model.spec}")

Model ID   : glm
Description: Standard Gaussian linear regression (GLM).
Y is generated as a linear combination of predictors X plus Gaussian noise.
Free parameters: regression weights (intercept + covariates).
Spec       : {'glm': {'linear': ['b0..bn']}}

In [2]:

# simulate data
nsubjects, nparams, ntrials = 80, 3, 100
true_params = np.random.randn(nsubjects, nparams)
X, Y = glm_model.sim(true_params, ntrials=ntrials)
all_data = [[X[i], Y[i]] for i in range(nsubjects)]

# simulate data nsubjects, nparams, ntrials = 80, 3, 100 true_params = np.random.randn(nsubjects, nparams) X, Y = glm_model.sim(true_params, ntrials=ntrials) all_data = [[X[i], Y[i]] for i in range(nsubjects)]

In [3]:

# fit and recover
model = EMModel(all_data=all_data, fit_func=glm_model.fit,
                param_names=[f"b{i}" for i in range(nparams)])
result = model.fit(verbose=0)

for param_idx, param_label in enumerate(range(nparams)):
    simulated_param = true_params[:,param_idx]
    estimated_param = model.outfit['params'][:,param_idx]
    ax = plotting.plot_scatter(simulated_param, f'Simulated b_{param_label}', 
                               estimated_param, f'Estimated b_{param_label}')

# fit and recover model = EMModel(all_data=all_data, fit_func=glm_model.fit, param_names=[f"b{i}" for i in range(nparams)]) result = model.fit(verbose=0) for param_idx, param_label in enumerate(range(nparams)): simulated_param = true_params[:,param_idx] estimated_param = model.outfit['params'][:,param_idx] ax = plotting.plot_scatter(simulated_param, f'Simulated b_{param_label}', estimated_param, f'Estimated b_{param_label}')

---

In [4]:

# try a logistic GLM
from pyem.models.glm import logit_model

print(f"Model ID   : {logit_model.id}")
print(f"Description: {logit_model.desc}")
print(f"Spec       : {logit_model.spec}")

# try a logistic GLM from pyem.models.glm import logit_model print(f"Model ID : {logit_model.id}") print(f"Description: {logit_model.desc}") print(f"Spec : {logit_model.spec}")

Model ID   : logit
Description: Standard logistic regression.
Y (0/1) is generated from a Bernoulli distribution with probability given by
the logistic (expit) link applied to a linear combination of predictors X.
Free parameters: regression weights (intercept + covariates).
Spec       : {'glm': {'linear': ['b0..bn']}, 'link': {'expit': []}}

In [5]:

# simulate data
nsubjects, nparams, ntrials = 80, 3, 100  # intercept + 2 predictors
rng = np.random.default_rng(0)
true_params = rng.normal(size=(nsubjects, nparams))

X, Y = logit_model.sim(true_params, ntrials=ntrials)
all_data = [[X[i], Y[i]] for i in range(nsubjects)]

# simulate data nsubjects, nparams, ntrials = 80, 3, 100 # intercept + 2 predictors rng = np.random.default_rng(0) true_params = rng.normal(size=(nsubjects, nparams)) X, Y = logit_model.sim(true_params, ntrials=ntrials) all_data = [[X[i], Y[i]] for i in range(nsubjects)]

In [6]:

# fit and recover
model = EMModel(
    all_data=all_data,
    fit_func=logit_model.fit,
    param_names=[f"b{i}" for i in range(nparams)],
)
result = model.fit(verbose=0)

# compare simulated vs estimated for each coefficient
for j in range(nparams):
    simulated_param = true_params[:, j]
    estimated_param = model.outfit["params"][:, j]
    ax = plotting.plot_scatter(
        simulated_param, f"Simulated b_{j}",
        estimated_param, f"Estimated b_{j}"
    )

# fit and recover model = EMModel( all_data=all_data, fit_func=logit_model.fit, param_names=[f"b{i}" for i in range(nparams)], ) result = model.fit(verbose=0) # compare simulated vs estimated for each coefficient for j in range(nparams): simulated_param = true_params[:, j] estimated_param = model.outfit["params"][:, j] ax = plotting.plot_scatter( simulated_param, f"Simulated b_{j}", estimated_param, f"Estimated b_{j}" )

---

In [7]:

# try an AR(1) GLM
from pyem.models.glm import glm_ar_model
from params import build_params

print(f"Model ID   : {glm_ar_model.id}")
print(f"Description: {glm_ar_model.desc}")
print(f"Spec       : {glm_ar_model.spec}")

# try an AR(1) GLM from pyem.models.glm import glm_ar_model from params import build_params print(f"Model ID : {glm_ar_model.id}") print(f"Description: {glm_ar_model.desc}") print(f"Spec : {glm_ar_model.spec}")

Model ID   : glm_ar
Description: Gaussian linear regression with an AR(1) autoregressive term
on the residuals: y_t = lin_t + phi * y_(t-1) + noise.
Free parameters: regression weights, phi (AR(1) coefficient, in (-1,1)).
Spec       : {'glm': {'linear': ['b0..bn'], 'ar1': ['phi']}}

In [8]:

# simulate data
nsubjects, nparams, ntrials = 80, 3, 100
true_params = np.random.randn(nsubjects, nparams - 1)  # regression weights (not registered - arbitrary length)

# phi (AR(1) coefficient) is registered, so draw it from the shared parameter registry (examples/params.py)
_, _, phi_col = build_params(["phi"], nsubjects)
true_params = np.hstack([true_params, phi_col])  # add phi
X, Y = glm_ar_model.sim(true_params, ntrials=ntrials)
all_data = [[X[i], Y[i]] for i in range(nsubjects)]

# simulate data nsubjects, nparams, ntrials = 80, 3, 100 true_params = np.random.randn(nsubjects, nparams - 1) # regression weights (not registered - arbitrary length) # phi (AR(1) coefficient) is registered, so draw it from the shared parameter registry (examples/params.py) _, _, phi_col = build_params(["phi"], nsubjects) true_params = np.hstack([true_params, phi_col]) # add phi X, Y = glm_ar_model.sim(true_params, ntrials=ntrials) all_data = [[X[i], Y[i]] for i in range(nsubjects)]

In [9]:

# fit and recover
param_names = [f"b{i}" for i in range(nparams-1)]+['phi']
model = EMModel(all_data=all_data, fit_func=glm_ar_model.fit,
                param_names=param_names,)
result = model.fit(verbose=0)

for param_idx, param_label in enumerate(param_names):
    simulated_param = true_params[:,param_idx]
    estimated_param = model.outfit['params'][:,param_idx]
    ax = plotting.plot_scatter(simulated_param, f'Simulated {param_label}', 
                               estimated_param, f'Estimated {param_label}')

# fit and recover param_names = [f"b{i}" for i in range(nparams-1)]+['phi'] model = EMModel(all_data=all_data, fit_func=glm_ar_model.fit, param_names=param_names,) result = model.fit(verbose=0) for param_idx, param_label in enumerate(param_names): simulated_param = true_params[:,param_idx] estimated_param = model.outfit['params'][:,param_idx] ax = plotting.plot_scatter(simulated_param, f'Simulated {param_label}', estimated_param, f'Estimated {param_label}')