import numpy as np from statsmodels.duration.hazard_regression import PHReg n = 300 # Sample size p = 5 # Number of covariates exog = np.random.normal(size=(n, p)) lin_pred = exog.sum(1) endog = -np.exp(-lin_pred)*np.log(np.random.uniform(size=n)) mod = PHReg(endog, exog) rslt = mod.fit() print rslt.summary() exog_rep = np.kron(exog, np.ones((4,1))) endog_rep = np.kron(endog, np.ones(4)) groups = np.kron(np.arange(n), np.ones(4)) mod_rep = PHReg(endog_rep, exog_rep) rslt_rep = mod_rep.fit() print rslt_rep.summary() mod_rep_a = PHReg(endog_rep, exog_rep) rslt_rep_a = mod_rep_a.fit(groups=groups) print rslt_rep_a.summary() n = 100 # Sample size p = 5 # Number of covariates exog = np.random.normal(size=(4*n, p)) lin_pred = exog.sum(1) endog = -np.exp(-lin_pred)*np.log(np.random.uniform(size=4*n)) groups = np.kron(np.arange(n), np.ones(4)) for method in "naive", "robust": if method == "robust": mod_rep_b = PHReg(endog, exog) rslt_rep_b = mod_rep_b.fit(groups=groups) else: mod_rep_b = PHReg(endog, exog) rslt_rep_b = mod_rep_b.fit() print rslt_rep_b.summary() print "\n" n = 300 # Sample size p = 5 # Number of covariates exog = np.random.normal(size=(n, p)) exog = np.kron(exog, np.ones((4,1))) lin_pred = exog.sum(1) endog = -np.exp(-lin_pred)*np.log(np.random.uniform(size=4*n)) groups = np.kron(np.arange(n), np.ones(4)) endog_a = endog.reshape((n,4)).mean(1) endog_a = np.kron(endog_a, np.ones(4)) # Adjust this parameter to control the amount of within-subject dependence # 0 = perfect dependence, 1 = no dependence a = 0.5 endog = a*endog + (1-a)*endog_a for method in "robust", "naive": if method == "robust": mod_rep_i = PHReg(endog, exog) rslt_rep_i = mod_rep_i.fit(groups=groups) else: mod_rep_i = PHReg(endog, exog) rslt_rep_i = mod_rep_i.fit() print rslt_rep_i.summary() print "\n"