Error r(430) when applying gammafit

[Tommy Wilkinson]

Hi there, I've installed gammafit successfully and generated the alpha and beta coefficients on expenditure data for my set of n=241 variables. When I attempt to use the gamma fit for a subset of the data (n=70 variables) using the "if" function I encountered an error that states "could not calculate numerical derivatives, missing values encountered r(430)".

To clarify: My variable of interest is "Total_costs". Within the dataset I have 3 age ranges (represented as binary variables) and I would like to generate the gamma distribution for each of the age ranges individually. Applying "gammafit Total_costs" generates the alpha and beta coefficients ok for the whole sample.

Applying "gammafit Total_costs if Age0_1 == 1" gives the error. I've looked at the help and seems like the "if" is allowed with the gammafit so I think I might be specifying it incorrectly? Assistance much appreciated.

[Nick Cox]

gammafit is from SSC. You are quite correct that if qualifiers (if is not a function) are allowed with gammafit. When you fit to a subset of 70 observations (not variables, presumably???) you are getting an error from the maximum likelihood engine which does almost all the work.

I suggest that you post that subset here as a data example. Then we can see whether a fit can be achieved.

[Nick Cox]

Here's something to try. A heuristic is that to a good approximation the cube root of a gamma-distributed variable is more nearly normally distributed.

So

gen curt_cost = Total_costs^(1/3)

qnorm curt_costs if Age0_1 == 1

A normal distribution would show on such a plot as essentially a linear configuration.

What do you see? I can't see your data but my wild guess is that you have one or both of two things:

1. A spike or cluster of values at or very near zero

2. Wild outliers that even the cube root can't tame.

It could be that such awkwardness are swamped in the total by the other 171 values.

[Tommy Wilkinson]

Many thanks for such a quick and useful response and apologies for my poor use of basic terminology!

For some further background context, my aim here is to generate cost distributions that I can use for a probabilistic sensitivity analysis in a cost effectiveness analysis of a health intervention. I've commonly seen gamma distributions applied for similar sensitivity analysis and my instinct is that the gamma distribution would be appropriate as the cost distribution will be non-zero with an expected long tail (i.e a few patients with very high costs).

See the data below. Note I also have 2 other age bands of Age2_4 and Age5_15 and the following code all returns the same error:
gammafit Total_costs if Age0_1== 1
gammafit Total_costs if Age2_4 == 1
gammafit Total_costs if Age5_15 == 1


I'm aiming to define the distribution for each age band because I would like to be able to run a PSA uniquely for each age band in my cost effectiveness analysis.

Total_cost | Age0_1
s | 0 1 | Total
-----------+----------------------+----------
840.94 | 1 0 | 1
873.79 | 1 0 | 1
875.15 | 1 0 | 1
919.47 | 1 0 | 1
951.73 | 1 0 | 1
958.40 | 1 0 | 1
959.36 | 1 0 | 1
960.16 | 1 0 | 1
986.95 | 1 0 | 1
987.78 | 1 0 | 1
996.68 | 1 0 | 1
1057.05 | 1 0 | 1
1062.21 | 1 0 | 1
1069.07 | 1 0 | 1
1081.00 | 1 0 | 1
1102.15 | 1 0 | 1
1115.66 | 1 0 | 1
1127.49 | 1 0 | 1
1133.69 | 0 1 | 1
1150.12 | 1 0 | 1
1182.95 | 1 0 | 1
1224.02 | 1 0 | 1
1273.55 | 1 0 | 1
1283.76 | 0 1 | 1
1285.36 | 1 0 | 1
1286.51 | 1 0 | 1
1331.91 | 1 0 | 1
1348.67 | 0 1 | 1
1351.47 | 1 0 | 1
1360.55 | 0 1 | 1
1370.86 | 1 0 | 1
1412.17 | 1 0 | 1
1459.92 | 1 0 | 1
1469.26 | 1 0 | 1
1471.67 | 1 0 | 1
1478.27 | 1 0 | 1
1500.75 | 0 1 | 1
1510.08 | 1 0 | 1
1515.71 | 1 0 | 1
1544.86 | 1 0 | 1
1545.24 | 1 0 | 1
1572.90 | 1 0 | 1
1574.83 | 1 0 | 1
1577.28 | 1 0 | 1
1589.85 | 1 0 | 1
1615.87 | 1 0 | 1
1635.76 | 1 0 | 1
1641.58 | 1 0 | 1
1686.70 | 1 0 | 1
1806.43 | 1 0 | 1
1808.56 | 0 1 | 1
1841.79 | 0 1 | 1
1852.90 | 0 1 | 1
1881.40 | 1 0 | 1
1892.53 | 1 0 | 1
1914.02 | 1 0 | 1
1966.99 | 0 1 | 1
1995.47 | 1 0 | 1
2053.90 | 1 0 | 1
2062.50 | 1 0 | 1
2107.16 | 1 0 | 1
2112.83 | 1 0 | 1
2139.33 | 1 0 | 1
2156.89 | 1 0 | 1
2161.72 | 0 1 | 1
2161.93 | 1 0 | 1
2169.77 | 1 0 | 1
2171.98 | 0 1 | 1
2185.00 | 0 1 | 1
2232.11 | 1 0 | 1
2330.66 | 0 1 | 1
2336.66 | 0 1 | 1
2369.32 | 0 1 | 1
2480.34 | 1 0 | 1
2507.15 | 0 1 | 1
2554.45 | 0 1 | 1
2566.09 | 1 0 | 1
2622.79 | 1 0 | 1
2792.47 | 0 1 | 1
2815.12 | 0 1 | 1
2851.25 | 0 1 | 1
2906.82 | 1 0 | 1
2997.16 | 1 0 | 1
3017.98 | 1 0 | 1
3044.80 | 0 1 | 1
3065.89 | 1 0 | 1
3162.83 | 1 0 | 1
3204.04 | 1 0 | 1
3240.88 | 1 0 | 1
3284.02 | 0 1 | 1
3301.45 | 0 1 | 1
3349.33 | 1 0 | 1
3399.42 | 1 0 | 1
3419.43 | 1 0 | 1
3447.68 | 1 0 | 1
3450.10 | 1 0 | 1
3521.58 | 1 0 | 1
3562.01 | 0 1 | 1
3596.46 | 1 0 | 1
3634.63 | 1 0 | 1
3668.84 | 0 1 | 1
3682.07 | 1 0 | 1
3698.84 | 1 0 | 1
3705.35 | 1 0 | 1
3719.95 | 1 0 | 1
3790.33 | 1 0 | 1
3854.42 | 1 0 | 1
3901.47 | 1 0 | 1
3916.04 | 0 1 | 1
3962.00 | 1 0 | 1
4016.85 | 1 0 | 1
4086.93 | 0 1 | 1
4217.28 | 1 0 | 1
4287.16 | 1 0 | 1
4308.35 | 1 0 | 1
4412.71 | 1 0 | 1
4485.77 | 1 0 | 1
4513.05 | 1 0 | 1
4712.39 | 1 0 | 1
4841.98 | 1 0 | 1
4888.24 | 0 1 | 1
5205.00 | 1 0 | 1
5233.51 | 0 1 | 1
5374.63 | 0 1 | 1
5376.61 | 0 1 | 1
5397.57 | 0 1 | 1
5518.36 | 1 0 | 1
5560.54 | 1 0 | 1
5631.12 | 1 0 | 1
5861.33 | 1 0 | 1
5921.43 | 1 0 | 1
6044.04 | 0 1 | 1
6099.11 | 0 1 | 1
6201.15 | 1 0 | 1
6290.74 | 1 0 | 1
6315.87 | 0 1 | 1
6608.59 | 1 0 | 1
6913.56 | 1 0 | 1
7066.16 | 0 1 | 1
7090.58 | 1 0 | 1
7348.97 | 1 0 | 1
7384.73 | 1 0 | 1
7509.43 | 1 0 | 1
7523.98 | 1 0 | 1
7659.87 | 1 0 | 1
7733.33 | 0 1 | 1
8596.18 | 0 1 | 1
8761.24 | 0 1 | 1
8933.79 | 0 1 | 1
9529.65 | 1 0 | 1
9716.22 | 1 0 | 1
9864.12 | 1 0 | 1
9999.93 | 1 0 | 1
10117.88 | 1 0 | 1
10126.89 | 1 0 | 1
10395.94 | 1 0 | 1
10527.01 | 1 0 | 1
10529.58 | 1 0 | 1
10746.43 | 1 0 | 1
11627.69 | 1 0 | 1
12094.39 | 0 1 | 1
12828.42 | 0 1 | 1
12965.15 | 1 0 | 1
13067.21 | 1 0 | 1
13518.18 | 1 0 | 1
13523.40 | 1 0 | 1
13787.47 | 0 1 | 1
13832.68 | 1 0 | 1
14044.52 | 1 0 | 1
14076.46 | 1 0 | 1
14215.62 | 0 1 | 1
14869.36 | 1 0 | 1
14873.04 | 1 0 | 1
15894.93 | 1 0 | 1
17154.21 | 0 1 | 1
17262.72 | 0 1 | 1
17267.08 | 0 1 | 1
17515.18 | 1 0 | 1
17820.52 | 0 1 | 1
18278.79 | 1 0 | 1
18639.33 | 1 0 | 1
18855.95 | 1 0 | 1
18877.68 | 1 0 | 1
19005.95 | 1 0 | 1
20091.61 | 1 0 | 1
20511.24 | 0 1 | 1
21009.26 | 0 1 | 1
21750.45 | 1 0 | 1
21859.29 | 1 0 | 1
22302.93 | 0 1 | 1
22386.16 | 1 0 | 1
22782.32 | 1 0 | 1
22865.15 | 0 1 | 1
23082.14 | 0 1 | 1
23343.00 | 0 1 | 1
23593.33 | 1 0 | 1
23843.12 | 1 0 | 1
24035.77 | 1 0 | 1
24208.26 | 0 1 | 1
25156.90 | 1 0 | 1
25831.11 | 0 1 | 1
26244.39 | 1 0 | 1
26781.97 | 0 1 | 1
26782.72 | 0 1 | 1
26800.01 | 1 0 | 1
26894.74 | 1 0 | 1
27008.16 | 1 0 | 1
27323.92 | 1 0 | 1
27580.36 | 1 0 | 1
28056.20 | 1 0 | 1
28553.72 | 1 0 | 1
28654.80 | 1 0 | 1
28849.63 | 1 0 | 1
28875.92 | 0 1 | 1
28879.64 | 1 0 | 1
29042.20 | 1 0 | 1
29084.09 | 1 0 | 1
29556.68 | 1 0 | 1
29751.22 | 0 1 | 1
30241.80 | 1 0 | 1
31138.93 | 1 0 | 1
32462.89 | 1 0 | 1
32818.35 | 0 1 | 1
32889.38 | 0 1 | 1
33463.65 | 1 0 | 1
34810.34 | 1 0 | 1
34963.24 | 0 1 | 1
35311.84 | 0 1 | 1
36004.36 | 0 1 | 1
36132.68 | 0 1 | 1
36295.77 | 1 0 | 1
36891.76 | 1 0 | 1
39587.00 | 1 0 | 1
42069.99 | 0 1 | 1
42610.43 | 1 0 | 1
42828.29 | 0 1 | 1
44377.64 | 0 1 | 1
45666.76 | 1 0 | 1
48905.06 | 0 1 | 1
60790.59 | 1 0 | 1
69829.53 | 1 0 | 1
-----------+----------------------+----------
Total | 171 70 | 241

I've also generated the graph as you suggested and this is what I got:


Note the Total_cost is the sum of all the other types of costs, and variable "Cost_Pharmaceuticals" is one type. When using code:

gammafit Cost_Pharmaceuticals if Age0_1 == 1

then it produces the coefficients successfully, so I imagine that you are correct that there is something about the clustering or outliers within "Total_costs".

Any insights much appreciated. Thanking you once again.

[Tommy Wilkinson]

And here is the graph noted above

[Tommy Wilkinson]

[Nick Cox]

Thanks for the data example. Please study FAQ Advice #12 https://www.statalist.org/forums/help#stata for

* Using dataex. Your listing required some engineering to be readable helpfully into Stata,

* Why posting .gph attachments doesn't work well. The forum software wasn't written by StataCorp and doesn't understand .gph file format.

I got gammafit to play by using technique(bhhh). I don't know more than I once tried it on a whim when facing difficulties like yours and it came up with results.

The resulting fit is a case of seen better, seen worse.

I used qgamma from SSC. for a quantile-quantile plot. Superimposing estimated density on a histogram would be interesting, but is more than i have time for right now.

* Example generated by -dataex-. For more info, type help dataex
clear
input double Cost byte Age0_1
  840.94 1
  873.79 1
  875.15 1
  919.47 1
  951.73 1
   958.4 1
  959.36 1
  960.16 1
  986.95 1
  987.78 1
  996.68 1
 1057.05 1
 1062.21 1
 1069.07 1
    1081 1
 1102.15 1
 1115.66 1
 1127.49 1
 1133.69 0
 1150.12 1
 1182.95 1
 1224.02 1
 1273.55 1
 1283.76 0
 1285.36 1
 1286.51 1
 1331.91 1
 1348.67 0
 1351.47 1
 1360.55 0
 1370.86 1
 1412.17 1
 1459.92 1
 1469.26 1
 1471.67 1
 1478.27 1
 1500.75 0
 1510.08 1
 1515.71 1
 1544.86 1
 1545.24 1
  1572.9 1
 1574.83 1
 1577.28 1
 1589.85 1
 1615.87 1
 1635.76 1
 1641.58 1
  1686.7 1
 1806.43 1
 1808.56 0
 1841.79 0
  1852.9 0
  1881.4 1
 1892.53 1
 1914.02 1
 1966.99 0
 1995.47 1
  2053.9 1
  2062.5 1
 2107.16 1
 2112.83 1
 2139.33 1
 2156.89 1
 2161.72 0
 2161.93 1
 2169.77 1
 2171.98 0
    2185 0
 2232.11 1
 2330.66 0
 2336.66 0
 2369.32 0
 2480.34 1
 2507.15 0
 2554.45 0
 2566.09 1
 2622.79 1
 2792.47 0
 2815.12 0
 2851.25 0
 2906.82 1
 2997.16 1
 3017.98 1
  3044.8 0
 3065.89 1
 3162.83 1
 3204.04 1
 3240.88 1
 3284.02 0
 3301.45 0
 3349.33 1
 3399.42 1
 3419.43 1
 3447.68 1
  3450.1 1
 3521.58 1
 3562.01 0
 3596.46 1
 3634.63 1
 3668.84 0
 3682.07 1
 3698.84 1
 3705.35 1
 3719.95 1
 3790.33 1
 3854.42 1
 3901.47 1
 3916.04 0
    3962 1
 4016.85 1
 4086.93 0
 4217.28 1
 4287.16 1
 4308.35 1
 4412.71 1
 4485.77 1
 4513.05 1
 4712.39 1
 4841.98 1
 4888.24 0
    5205 1
 5233.51 0
 5374.63 0
 5376.61 0
 5397.57 0
 5518.36 1
 5560.54 1
 5631.12 1
 5861.33 1
 5921.43 1
 6044.04 0
 6099.11 0
 6201.15 1
 6290.74 1
 6315.87 0
 6608.59 1
 6913.56 1
 7066.16 0
 7090.58 1
 7348.97 1
 7384.73 1
 7509.43 1
 7523.98 1
 7659.87 1
 7733.33 0
 8596.18 0
 8761.24 0
 8933.79 0
 9529.65 1
 9716.22 1
 9864.12 1
 9999.93 1
10117.88 1
10126.89 1
10395.94 1
10527.01 1
10529.58 1
10746.43 1
11627.69 1
12094.39 0
12828.42 0
12965.15 1
13067.21 1
13518.18 1
 13523.4 1
13787.47 0
13832.68 1
14044.52 1
14076.46 1
14215.62 0
14869.36 1
14873.04 1
15894.93 1
17154.21 0
17262.72 0
17267.08 0
17515.18 1
17820.52 0
18278.79 1
18639.33 1
18855.95 1
18877.68 1
19005.95 1
20091.61 1
20511.24 0
21009.26 0
21750.45 1
21859.29 1
22302.93 0
22386.16 1
22782.32 1
22865.15 0
23082.14 0
   23343 0
23593.33 1
23843.12 1
24035.77 1
24208.26 0
 25156.9 1
25831.11 0
26244.39 1
26781.97 0
26782.72 0
26800.01 1
26894.74 1
27008.16 1
27323.92 1
27580.36 1
 28056.2 1
28553.72 1
 28654.8 1
28849.63 1
28875.92 0
28879.64 1
 29042.2 1
29084.09 1
29556.68 1
29751.22 0
 30241.8 1
31138.93 1
32462.89 1
32818.35 0
32889.38 0
33463.65 1
34810.34 1
34963.24 0
35311.84 0
36004.36 0
36132.68 0
36295.77 1
36891.76 1
   39587 1
42069.99 0
42610.43 1
42828.29 0
44377.64 0
45666.76 1
48905.06 0
60790.59 1
69829.53 1
end


gammafit Cost if Age0_1 == 1 ,  technique(bhhh) 

qgamma Cost if Age0_1 == 1, param(.8569521 12095.76)

[Tommy Wilkinson]

Thanks a million, absolute lifesaver. Code working well and generating results with al 3 age brackets. Out of interest here is the histogram using code:
histogram Total_costs if Age0_1==1, width(400) addplot(function gammaden(e(alpha), e(beta), 0, x), ra(0 40000))

[Nick Cox]

Thanks for the report. I suggest a bigger bin width on your histogram or comparison with a kernel density estimate.

[Tommy Wilkinson]

Thanks again Nick, very much appreciated. Yes I've widened bin width on the histogram to 1000 and much improved. See also the kernel density plot as suggested, which does seem to support the gamma distribution

[Tommy Wilkinson]

Plot as noted

[Tommy Wilkinson]

One further question if that's OK? The calculation of the mean from the gammafit is straightforward as = alpha*beta coefficients. I have the standard error and 95% confidence interval for the alpha and beta coefficients, but can you advise code for how I would calculate a credible interval for the mean?

[Nick Cox]

"Support" is a stronger word than I would use here. It depends on dismissing the secondary mode around 25000 as a sampling quirk. Also, IIUC #8 implies a mode at zero but the density estimate in #12 doesn't.

There are several ways to get a confidence interval for the mean and I don't have an informed take on any being best, What's a credible interval if it's not a confidence interval? You could also try a generalised linear model with a gamma family,

[Tommy Wilkinson]

Thanks again for this Nick I really appreciate it. I'm working my way through applying a generalised linear model with a gamma family as you suggest above. I've run the GLM model of total costs against one continuous variable treatment length (in months) using code and output:

. glm Total_costs Trmt_duration_month, family(gamma) link(log) vce(bootstrap)

Note also the output from the gammafit on Total_costs:

Apologies this feels like a very basic question, but how would I then calculate (i.e what would be the code) for a confidence interval around the mean of Total_costs from this output? i.e as I know the distribution is not normally distributed I assume it would not be appropriate to use:
. ci means Total_costs

And I see from the Stata help file on ci there is alternative syntax for confidence interval for means for a Poisson distribution, but I couldn't find a comparable approach for confidence interval for means for a gamma distribution.

Thanks again.

[Nick Cox]

There are two basic tasks being conflated here. Fitting a gamma distribution on its own; and fitting a model with predictors. A good result for one doesn't rule out a poor result for the other, whichever way you take it.

You could be model or distribution agnostic and get some of bootstrap confidence interval for your mean.

The following sequence shows some of the possibilities.

sysuse auto, clear 

ci mean price

ci mean price, poisson

glm price , f(gamma) vce(robust) link(identity)

glm price , f(poisson) vce(robust) link(identity)

bootstrap r(mean), r(1000) nodots : su price, meanonly

estat bootstrap, all