yaesoubilab / simpy Goto Github PK

functions get_CI and get_PI in ICER_indp

Note that all assert statements will be removed when you run the code. They are only active if run your code in the debug mode. So instead of assert, it's better to check the necessary condition and if it's not satisfied use raise statement.

Hi Zongbo,

I pushed the changed I made. Make sure to get the latest version of the repository and switch to the "TestingOptimization" branch to see these changes. You will notice that now that we are using common random number seeds to calculate derivatives, the convergence has improved substantially.

I've also added two separate files to test both constrained and unconstrained scenarios.

We should use the test functions to make sure that everything is working correctly. So I wonder if you could do the following:

1. Write a function similar to plot_x_irs() (in Optimization.py file) to visualize how the derivative is changing over iterations. Note that derivatives are stored in self.itr_nDf.
2. Write a function to create single plot to visualize how the harmonic step size and the derivative step size are changing over iterations. Note that value of these step sizes are stored in step_df and step_move for each iteration.

And let me know whenever you push your changes and I'll review them.

In function build_CE_table (in file EconEvalClasses.py), there is variable named output_estimates and a class member with the same name self.output_estimates. I don't think we are using self.output_estimates anywhere outside this function. So it is better to convert it to a regular variable. I was going to change it but didn't want to mess it up! I think output_estimates and self.output_estimates are different.

Could you add comments to explain better how the maximum likelihood algorithm works for gamma Poisson?

Whenever you get a chance, would you please write the code (and the test code) to fit these distributions to some synthesized data?

1. Exponential
2. Beta (note that the data may not necessarily fall within 0 and 1, so we should take min and max as optional arguments for fit_beta(). If they are not provided, they the code needs to find the min and max from the data.
3. BetaBinomial
4. Binomial
5. Empirical (I guess for this, we just need to return the frequency of each observation)
6. Gamma
7. GammaPoisson
8. Geometric
9. JohnsonSb
10. JohnsonSu
11. LogNormal
12. NegativeBinomial
13. Normal
14. Triangular
15. Uniform
16. UniformDiscrete
17. Weibull

Some must be easy some more challenging. So let me know if you have any questions!

In the file TestEconEval_CEA.py, in line 42, I get an error if instead of
print(myCEA.build_CE_table(ce.Interval.PREDICTION))
I use
print(myCEA.build_CE_table()).

In other words, the code doesn't like it when we ask for the CE table without intervals.

I imagine the red curve should be identical to the histogram. Do you agree? I think I know what is happening but let's discuss the best way to present this.

shouldn't it return n as well?

Could you add comments to provide more details about the MLE algorithm for negative binomial?

This is mainly a question. I thought AIC is always positive (but I am not sure at all). If that's true, then we are getting some negative AICs.

In CETable, would it be possible to add the 1000 separator (",") to all the numbers if necessary?

function show_CE_plane

I have added three new classes at the end of StatisticalClasses.py to make inference about the random variable (X-Y)/X.

Whenever you have a chance, would you please implement these classes (they will look similar to _RatioStat , RatioStatPaired, and RatioStatIndp classes).

And note that since (X-Y)/X = 1-Y/X, you can tie these new classes to the corresponding ratio classes.

Use "" instead of '

Could you add comments to explain in more details how the maximum likelihood estimates for beta binomial are obtained? Function fit_beta_binomial in FittingProbDist_ML.py.

And estimate n from data.

I know you asked me before and I answered no but now I think it's often needed to get the shifted strategies that are on the frontier and not on the frontier. So would you please implement the following two methods in the CEA class?

def get_shifted_strategies_on_frontier(self):
def get_shifted_strategies_not_on_frontier(self):

And remember both need to return a list (not data frame).

Raise error if not enough observations are provided to calculate the confidence interval.

Could show_CE_plane take two optional arguments (x_range=None, y_range=None) to take the range of x and y axes? You can assume that x_range and y_range are provided as list of [l, u].

I tried to do it myself but then the axis ticks didn't show up correctly.

1. To initialize CEA class, strategies are provided in list format. Could we modify the code so that self._strategiesOnFrontier is also the list of strategies that are on the frontier?
2. And could we also assign self._strategiesNotOnFrontier to the list of strategies that are not on the frontier?

Would it be possible to format the numbers on the y-axis of the cost-effectiveness figure (in EconEvalClasses.py) such that they are rounded to an integer and have 1000 separator (,) if needed?

When using the method of moments (FittingProbDist_MM.py), shouldn't get_binomial_parameters return n as well?

[Make sure to pull the latest version]

The bootstrap algorithm for DifferenceStatIndp and RatioStatIndp is different from the bootstrap algorithm for SummaryStat.

Also, when testing on data, the bootstrap CI and t-CI in SummaryStat yields almost the same number (as expected) that gives me the feeling that the bootstrap version in SummarStat is correct.

What do you think?

First, notice that I changed the set up in TestCEPlot.py so that the strategy s0 is not on the origin any more. This is more realistic. So now, let's assume that the first strategy in the list strategies always represents the "base" strategy, and so we need to move the origin to this strategy:

The the __init__ function of class CEA should be modified such that it takes a new argument if_paired.

if_paired = False then subtract E[cost] and E[effect] of the first strategy from the E[cost] and E[effect] of all strategies so that the first strategy is now on the origin.
And
if_paired = False the subtract cost_i and effect_i of the first strategy from cost_i and effect_i of all strategies so that the first strategy is now on the origin.

using non-default values

This is hopefully easy. Would it be possible to revise the show_CE_plane method in EconEvalClasses.py file so that when the range of x or y-axis is less than 10, we let python to automatically format the numbers and ticks (so that the decimal points show up as needed)?

I had a student whose cost and effect estimates were less than 1, and when she used the show_CE_plane method, all tick labels were 0.

Thanks so much.

Modify so that it also generate sample path for non-absorbing health states.

For independent variable x, y, the ratio variable have:
E(x/y) = E(x)/E(y)
var(x/y) = E(x^2)/E(y^2) - E(x)^2/E(y)^2

But the distribution of E(x/y) is different from the distribution of E(x)/E(y)
(similar reason applied to E(x-y) and E(x)-E(y))

if both have mu=0, x/y has a Cauchy distribution
https://en.wikipedia.org/wiki/Ratio_distribution#Gaussian_ratio_distribution

For the class CEA in EconEvalClasses.py, could you add "Frontier" to the legend of the CE figure when the show_legend argument of show_CE_plane is True?

Decision and Chance nodes should not accept utility. Only Terminal nodes can have utility.

Just to confirm, for the empirical distribution, is the input data like:
(1, 1, 2, 3, 2)
or
(1, 40%), (2, 40%), (3, 20%) ?
Thanks a lot!

I can't find useful information regarding the Johnson SI distribution. Do you mind give me some references about it? Thank you!

for when Paired and for interval == Interval.CONFIDENCE

[Pull first]

In the cost-effectiveness table, the columns E[Cost] and E[Effect] should show the average cost of each strategies but it current shows the average incremental cost and effect with respect to the first strategy.

Would you include some sample code to fit probability distributions to some randomly generated data? I have added ProbDistFitting.py for this purpose. The sample code should show the histogram of the data on which the fitted probability distribution is overlaid.

Let's start with Poisson, and Gamma for now.

1. get from the user number of significant digits to round the Cost, Health and ICER to.
2. print table to the console
3. print table to a csv file

I think this looks great. You did an excellent job.
We just need to make some changes to make it more readable and easy to use:

1. The class Strategy represents each of the point to be shown on the CE figure. To initialize a Strategy, we also need to give it a name. So please modify the __init__ function to: def __init__(self, name, cost, effect), where name is the name of this strategy, cost is the expected cost of this strategy, and effect is the expected effect of this strategy.

2. Therefore, class Strategy should have these members: self.name, self.cost, self.effect, and self.ifDominated.

3. Now to initialize CEA class, we should take the list strategies which contains all strategies.

4. Please modify FindFrontier method so that it returns the list of strategies that will be on the frontier (do not show the figure yet). And note that this method should also update the list self.strategies so that you can use this later one to generate the figure and the table.

5. Add a new method ShowCEPlane that displays the cost-effectiveness figure.

6. On the figure, please make the x and y axes a little thinner, and the frontier Black.

In SamplePathClass.py, could you please make the following changes to the function graph_sample_paths?

1. Use the argument transparency to change the transparency of plotted curves.
2. If the argument if_same_color is set to True, paint the sample paths the same color.

The tests for this class are located at TestSamplePath.py.