Forecasts from LGBExitModeler().forecast() do not sum to 1 within Individual-time period. However, they are close. A manual normalization reduces the error by half, as shown in code. However, a manual normalization could obscure other bugs and make debugging harder in the future.
from fife.lgb_modelers import LGBSurvivalModeler
from fife.lgb_modelers import LGBExitModeler
from fife.processors import PanelDataProcessor
from fife.utils import make_results_reproducible
from pandas import concat, date_range, read_csv, to_datetime, Index, DataFrame
#import pandas as pd
#import numpy as np
from numpy import array, repeat, full, isclose
# import warnings
# warnings.filterwarnings('ignore')
SEED = 9999
make_results_reproducible(SEED)
data = read_csv("https://www.dl.dropboxusercontent.com/s/3tdswu2jfgwp4xw/REIGN_2020_7.csv?dl=0")
data.head()
data["country-leader"] = data["country"] + ": " + data["leader"]
data["year-month"] = data["year"].astype(int).astype(str) + data["month"].astype(int).astype(str).str.zfill(2)
data["year-month"] = to_datetime(data["year-month"], format="%Y%m")
data = concat([data[["country-leader", "year-month"]],
data.drop(["ccode", "country-leader", "leader", "year-month"],
axis=1)],
axis=1)
data.head()
total_obs = len(data)
data = data.drop_duplicates(["country-leader", "year-month"], keep="first")
# n_duplicates = total_obs - len(data)
# print(f"{n_duplicates} observations with a duplicated identifier pair deleted.")
data_processor = PanelDataProcessor(data=data)
data_processor.build_processed_data()
# data_processor.data.head()
mdata = data_processor.data.copy()
# Assign unique ID based off broken panels
gaps = mdata.groupby("country-leader")["_period"].shift() < mdata["_period"] - 1
spells = gaps.groupby(mdata["country-leader"]).cumsum()
mdata['country-leader-spell'] = mdata['country-leader'] + array(spells, dtype = 'str')
## Create exit status for dissolution, change of government type, and same government type
mdata['outcome'] = 'no_exit'
for i in mdata['country-leader-spell'].unique():
# Boolean index for the unique country-leader-spell observations
dex1 = mdata['country-leader-spell'] == i
# Boolean index for the last observation for the country-leader-spell observations
dex2 = dex1 & (mdata['year-month'] == mdata.loc[dex1,'year-month'].max())
# Boolean index for the last observation for the country-leader-spell observations unless it is censored
dex3 = dex2 & (mdata['year-month'] != mdata['year-month'].max())
# Only assign exit status if not censored
if dex3.sum() > 0:
futuredex = (mdata['year-month'] >= mdata.loc[dex3,'year-month'].values[0]) & (mdata['country'] == mdata.loc[dex3,"country"].values[0]) & ~dex1
if futuredex.sum() == 0:
mdata.loc[dex2,'outcome'] = "Dissolution" # could alternatively use dex1
else:
oldgovernment = mdata.loc[dex3,'government'].values[0]
futuremdata = mdata[futuredex]
newgovernment = futuremdata.loc[futuremdata['year-month'] == futuremdata['year-month'].min(),'government'].values[0]
if oldgovernment == newgovernment:
mdata.loc[dex2,'outcome'] = "Same_government_type" # could alternatively use dex1
else:
mdata.loc[dex2,'outcome'] = "Change_of_government_type" # could alternatively use dex1
# We no longer need country-leader-spell
mdata = mdata.drop('country-leader-spell', axis = 1)
# The outcome variable must be a category type
mdata['outcome'] = mdata['outcome'].astype('category')
# Obtain probability of type of exit conditional on exit
exit_modeler = LGBExitModeler(data=mdata, exit_col = 'outcome')
exit_modeler.build_model(parallelize=False)
exit_modeler_forecasts = exit_modeler.forecast()
# Do the outcomes always sum to 1?
summary1a = DataFrame(index = exit_modeler_forecasts.index.unique(), columns = exit_modeler_forecasts.columns[:-1])
for i in exit_modeler_forecasts.index.unique():
summary1a.loc[i,:] = exit_modeler_forecasts.loc[i, exit_modeler_forecasts.columns != 'Future outcome'].sum(axis = 0)
(summary1a!=1).sum().sum()/(summary1a.shape[0]*summary1a.shape[1]) # Should be 0
abs(1-summary1a).sum().sum()/(summary1a.shape[0]*summary1a.shape[1]) # should be 0
# However, they are "close"
tmp = summary1a.to_numpy(dtype='float64')
isclose(tmp,full(tmp.shape,1)).sum()/(tmp.shape[0] * tmp.shape[1])
# Slightly improved normalization
exit_modeler_forecasts_normalized = exit_modeler_forecasts.copy()
exit_modeler_forecasts_normalized.iloc[:,:-1] = exit_modeler_forecasts_normalized.iloc[:,:-1].values / DataFrame(repeat(summary1a.values, repeats = mdata['outcome'].nunique(), axis = 0)).values
# Do the outcomes always sum to 1?
summary1a_normalized = DataFrame(index = exit_modeler_forecasts_normalized.index.unique(), columns = exit_modeler_forecasts_normalized.columns[:-1])
for i in exit_modeler_forecasts.index.unique():
summary1a_normalized.loc[i,:] = exit_modeler_forecasts_normalized.loc[i, exit_modeler_forecasts_normalized.columns != 'Future outcome'].sum(axis = 0)
(summary1a_normalized!=1).sum().sum()/(summary1a_normalized.shape[0]*summary1a_normalized.shape[1]) # Should be 0
abs(1-summary1a_normalized).sum().sum()/(summary1a_normalized.shape[0]*summary1a_normalized.shape[1]) # should be 0
# half as many different than one and absolute error is cut in half