Я думаю, входная переменная x для group_rand_coef_data неизвестен при запуске этой модели смешанных эффектов в этой строке gpb.GPModel(group_data=group_data, group_rand_coef_data=x, в реальной жизни, в отличие от этого кода руководства (ниже), где x был известен/ создано человеком.
В реальной жизни, когда x неизвестен, что я могу использовать для group_rand_coef_data?
Код: Выделить всё
import gpboost as gpb
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
plt.style.use('ggplot')
def simulate_response_variable(lp, rand_eff, likelihood):
"""Function that simulates response variable for various likelihoods"""
n = len(rand_eff)
if likelihood == "gaussian":
xi = 0.1**0.5 * np.random.normal(size=n) # error term, variance = 0.1
y = lp + rand_eff + xi
return y
likelihood = "gaussian"
"""
Grouped random effects
"""
# --------------------Simulate data----------------
# Single-level grouped random effects
n = 1000 # number of samples
m = 200 # number of categories / levels for grouping variable
group = np.arange(n) # grouping variable
for i in range(m):
group[int(i * n / m):int((i + 1) * n / m)] = i
np.random.seed(1)
b = 0.25**0.5 * np.random.normal(size=m) # simulate random effects, variance = 0.25
# Simulate linear regression fixed effects
X = np.column_stack((np.ones(n), np.random.uniform(size=n) - 0.5)) # design matrix / covariate data for fixed effect
beta = np.array([0, 2]) # regression coefficents
lp = X.dot(beta)
# Crossed grouped random effects and random slopes
group_crossed = group[np.random.permutation(n)-1] # grouping variable for crossed random effects
b_crossed = 0.25**0.5 * np.random.normal(size=m) # simulate crossed random effects
b_random_slope = 0.25**0.5 * np.random.normal(size=m)
x = np.random.uniform(size=n) # covariate data for random slope
rand_eff = b[group] + b_crossed[group_crossed] + x * b_random_slope[group]
rand_eff = rand_eff - np.mean(rand_eff)
y_crossed_random_slope = simulate_response_variable(lp=lp, rand_eff=rand_eff, likelihood=likelihood)
# --------------------Two crossed random effects and random slopes----------------
# Define and train model
group_data = np.column_stack((group, group_crossed))
gp_model = gpb.GPModel(group_data=group_data, group_rand_coef_data=x,
ind_effect_group_rand_coef=[1], likelihood=likelihood)
# 'ind_effect_group_rand_coef=[1]' indicates that the random slope is for the first random effect
gp_model.fit(y=y_crossed_random_slope, X=X, params={"std_dev": True})
gp_model.summary()
# Prediction
pred = gp_model.predict(group_data_pred=group_data, group_rand_coef_data_pred=x, X_pred=X)
# Obtain predicted (="estimated") random effects for the training data
all_training_data_random_effects = gp_model.predict_training_data_random_effects()
first_occurences_1 = [np.where(group==i)[0][0] for i in np.unique(group)]
pred_random_effects = all_training_data_random_effects.iloc[first_occurences_1,0]
pred_random_slopes = all_training_data_random_effects.iloc[first_occurences_1,2]
# Compare true and predicted random effects
plt.scatter(b, pred_random_effects, label="Random effects")
plt.scatter(b_random_slope, pred_random_slopes, label="Random slopes")
plt.legend()
plt.title("Comparison of true and predicted random effects")
plt.show(block=False)
Код: Выделить всё
=====================================================
Model summary:
Log-lik AIC BIC
-823.6 1659.2 1688.64
Nb. observations: 1000
Nb. groups: 200 (Group_1), 200 (Group_2)
-----------------------------------------------------
Covariance parameters (random effects):
Param. Std. dev.
Error_term 0.1044 0.0067
Group_1 0.1918 0.0263
Group_2 0.2709 0.0301
Group_1_rand_coef_nb_1 0.2288 0.0527
-----------------------------------------------------
Linear regression coefficients (fixed effects):
Param. Std. dev. z value P(>|z|)
Covariate_1 -0.0051 0.051 -0.0995 0.9208
Covariate_2 1.9975 0.048 41.6060 0.0000
=====================================================
Подробнее здесь: https://stackoverflow.com/questions/784 ... ts-gpmodel
