\(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\W}[1]{ \; #1 \; }\)
user_predict_mulcov.py#
View page sourcePredict Covariate Multiplier Values#
Purpose#
This examples demonstrates covariate multiplier predictions.
Problem Parameters#
The following values are used to simulate the data and set the priors for fitting the data:
age_list = [ 0.0, 100.0 ]
time_list = [ 1990.0, 2015.0 ]
def alpha_true(age, time) :
age_factor = 1.0 + (age - age_list[0]) / (age_list[1] - age_list[0])
time_factor = 1.0 + (time - time_list[0]) / (time_list[1] - time_list[0])
return 0.1 * age_factor * time_factor
iota_reference = 1e-2
income_reference = 1.0
Age and Time Values#
The reference value for iota is constant in age and time, but the value of the covariate multiplier that affects iota changes with age and time.
Rate Table#
The rate_table specifies that the only rate variable is iota for north_america. In addition, it specifies the smoothing for this rate has one grid point.
Mulcov Table#
The mulcov_table specifies that
the covariate alpha is a bilinear function of age and time.
In fact, it is equal to the function alpha_true defined
as one of the problem parameters.
Variables#
There are five model variables in this example:
iota_reference |
There is one variable corresponding to the reference value for iota(a,t) in north_america. |
alpha |
There are four variables corresponding to the rate_value covariate multiplier that affect iota . |
Integrand Table#
The integrand_table for this example includes
Sincidence and
mulcov .
The mulcov_0 integrand corresponds the value of alpha .
Data Table#
There are four measurements of Sincidence in the data_table , one for (age, time) pair corresponding to an alpha model variable. No noise is added to the measurements, and the prior on iota constrains it to the iota_reference .
Avgint Table#
There are four predictions of alpha requested by the avgint_table , one for each alpha model variable. The predictions are compared with the truth to see that the fit is prefect (there is no noise in the data).
Source Code#
# begin problem parameters
age_list = [ 0.0, 100.0 ]
time_list = [ 1990.0, 2015.0 ]
def alpha_true(age, time) :
age_factor = 1.0 + (age - age_list[0]) / (age_list[1] - age_list[0])
time_factor = 1.0 + (time - time_list[0]) / (time_list[1] - time_list[0])
return 0.1 * age_factor * time_factor
iota_reference = 1e-2
income_reference = 1.0
# end problem parameters
# ---------------------------------------------------------------------------
import sys
import os
import copy
import math
test_program = 'example/user/predict_mulcov.py'
if sys.argv[0] != test_program or len(sys.argv) != 1 :
usage = 'python3 ' + test_program + '\n'
usage += 'where python3 is the python 3 program on your system\n'
usage += 'and working directory is the dismod_at distribution directory\n'
sys.exit(usage)
print(test_program)
#
# import dismod_at
local_dir = os.getcwd() + '/python'
if( os.path.isdir( local_dir + '/dismod_at' ) ) :
sys.path.insert(0, local_dir)
import dismod_at
#
# change into the build/example/user directory
if not os.path.exists('build/example/user') :
os.makedirs('build/example/user')
os.chdir('build/example/user')
# ---------------------------------------------------------------------------
def example_db (file_name) :
def fun_income(a, t) :
return ('prior_income_value', 'prior_income_dage', 'prior_income_dtime')
def fun_iota_parent(a, t) :
return ('prior_iota_parent', None, None)
# ----------------------------------------------------------------------
# age table
#
# time table
node_table = [
{ 'name':'north_america', 'parent':'' },
]
#
# weight table:
weight_table = list()
# integrand table
integrand_table = [
{ 'name':'Sincidence' },
{ 'name':'mulcov_0' }
]
#
# covariate table:
covariate_table = [{ 'name':'income', 'reference': income_reference }]
#
# mulcov table: one multiplier with mulcov_id = 0
mulcov_table = [{
'covariate': 'income',
'type': 'rate_value',
'effected': 'iota',
'group': 'world',
'smooth': 'smooth_income'
} ]
#
# nslist_dict:
nslist_dict = dict()
# ----------------------------------------------------------------------
# prior_table
prior_table = [
{ # prior_iota_parent
'name': 'prior_iota_parent',
'density': 'uniform',
'lower': iota_reference,
'upper': iota_reference,
'mean': iota_reference,
},{ # prior_income_value
'name': 'prior_income_value',
'density': 'uniform',
'mean': 0.0,
'lower': -1.0,
'upper': +1.0,
},{ # prior_income_dage
'name': 'prior_income_dage',
'density': 'uniform',
'mean': 0.0,
},{ # prior_income_dtime
'name': 'prior_income_dtime',
'density': 'uniform',
'mean': 0.0,
}
]
# ----------------------------------------------------------------------
# smooth table
smooth_table = [
{ # smooth_iota_parent
'name': 'smooth_iota_parent',
'age_id': [ 0 ],
'time_id': [ 0 ],
'fun': fun_iota_parent
}, { # smooth_income
'name': 'smooth_income',
'age_id': [ 0 , 1],
'time_id': [ 0 , 1],
'fun': fun_income
}
]
# ----------------------------------------------------------------------
# rate table
rate_table = [{
'name': 'iota',
'parent_smooth': 'smooth_iota_parent',
} ]
# --------------------------------------------------------------------
# option_table
option_table = [
{ 'name':'parent_node_name', 'value':'north_america' },
{ 'name':'rate_case', 'value':'iota_pos_rho_zero' },
{ 'name':'quasi_fixed', 'value':'true' },
{ 'name':'max_num_iter_fixed', 'value':'30' },
{ 'name':'print_level_fixed', 'value':'0' },
{ 'name':'tolerance_fixed', 'value':'1e-10' },
]
# ----------------------------------------------------------------------
# data table:
data_table = list()
row = {
'node': 'north_america',
'subgroup': 'world',
'density': 'gaussian',
'weight': '',
'hold_out': False,
'integrand': 'Sincidence',
'income': 2.0,
}
for age in age_list :
for time in time_list :
alpha = alpha_true(age, time)
effect = alpha * ( row['income'] - income_reference )
iota = iota_reference * math.exp( effect )
row['age_lower'] = age
row['age_upper'] = age
row['time_lower'] = time
row['time_upper'] = time
row['meas_value'] = iota
row['meas_std'] = iota / 10.0
data_table.append( copy.copy(row) )
# ----------------------------------------------------------------------
# avgint table:
avgint_table = list()
# values that are the same for all rows
does_not_matter = 5.0
row = {
'node': 'north_america',
'subgroup': 'world',
'integrand': 'mulcov_0',
'weight': '',
'income': does_not_matter,
}
for age in age_list :
for time in time_list :
row['age_lower'] = age
row['age_upper'] = age
row['time_lower'] = time
row['time_upper'] = time
avgint_table.append( copy.copy(row) )
# ----------------------------------------------------------------------
# subgroup_table
subgroup_table = [ { 'subgroup':'world', 'group':'world' } ]
# ----------------------------------------------------------------------
# create database
dismod_at.create_database(
file_name,
age_list,
time_list,
integrand_table,
node_table,
subgroup_table,
weight_table,
covariate_table,
avgint_table,
data_table,
prior_table,
smooth_table,
nslist_dict,
rate_table,
mulcov_table,
option_table
)
# ===========================================================================
file_name = 'example.db'
example_db(file_name)
#
program = '../../devel/dismod_at'
dismod_at.system_command_prc([ program, file_name, 'init' ])
dismod_at.system_command_prc([ program, file_name, 'fit', 'fixed' ])
dismod_at.system_command_prc([ program, file_name, 'predict', 'fit_var' ])
# -----------------------------------------------------------------------
# connect to database
connection = dismod_at.create_connection(
file_name, new = False, readonly = True
)
predict_table = dismod_at.get_table_dict(connection, 'predict')
avgint_table = dismod_at.get_table_dict(connection, 'avgint')
node_table = dismod_at.get_table_dict(connection, 'node')
#
for row in predict_table :
avg_integrand = row['avg_integrand']
avgint_id = row['avgint_id']
age = avgint_table[avgint_id]['age_lower']
time = avgint_table[avgint_id]['time_lower']
truth = alpha_true(age, time)
assert abs( 1.0 - avg_integrand / truth ) < 1e-6
# -----------------------------------------------------------------------------
print('predict_mulcov.py: OK')