non-central t, power, R: [🇷 for BE/BA]
Dear all,
In a previous thread I asked what the power in a 2,2,2-BE trial is at N=6, CV=65% and T/R=95%.
The things is, Fartssie comes up with a negative reply. HS's implementation in R of the equations given by Julious (Statistics in Medicine, 2004 vol. 23, 1921-1986) also evaluates to a negative figure.
A central role here is played by the non-central t distribution. My guess was that we were torturing the algorithm, and I think this was HS's point as well (quote "I think we are pushing software to the limits, i.e., running into troubles of getting a reasonable value of the noncentral t-distribution (...)").
This I wanted to look a little into using R as my reference. Wikipedia tells us what the distribution's density function looks like. I made an independent implementation of this expression, plus an implementation of the integrated density function from t and downwards, which is basically the same as the function pt in R which is based on ASA243 according to the documentation.
The values from pt can be easily reproduced, it turns out. But I cannot "improve" them, although I have all the time in the world to integrate with ridiculously small stepsizes and tweak and fiddle. As it turns out in the concrete example, we are asking for probabilities from the noncentral t dist in a region of t-values plusminus 2.13 and noncentrality parameters 0.50 and -0.80. These are not at all extreme values. I think it sounds reasonable that "prob1" and "prob2" (using HS's terminology) are around 0.89 and 0.17 [try and plot the densities as function of t, eyeball the relative amount of curve left of the indicated values].
For this reason I am becoming a sceptic now. Either I am not able to make a better implementation than ASA243, or ...please forgive me... the equations given by Julious (and others) are not exact in certain situations (for example: I only see the issue when t2 < nc2), or something entirely else is in play.
Can anyone shed light on this?
Let me add: If the problem lies in ASA243 then I'd love to write one that gets it right, of course.)
Many thanks for any input and best regards,
EM.
In a previous thread I asked what the power in a 2,2,2-BE trial is at N=6, CV=65% and T/R=95%.
The things is, Fartssie comes up with a negative reply. HS's implementation in R of the equations given by Julious (Statistics in Medicine, 2004 vol. 23, 1921-1986) also evaluates to a negative figure.
A central role here is played by the non-central t distribution. My guess was that we were torturing the algorithm, and I think this was HS's point as well (quote "I think we are pushing software to the limits, i.e., running into troubles of getting a reasonable value of the noncentral t-distribution (...)").
This I wanted to look a little into using R as my reference. Wikipedia tells us what the distribution's density function looks like. I made an independent implementation of this expression, plus an implementation of the integrated density function from t and downwards, which is basically the same as the function pt in R which is based on ASA243 according to the documentation.
The values from pt can be easily reproduced, it turns out. But I cannot "improve" them, although I have all the time in the world to integrate with ridiculously small stepsizes and tweak and fiddle. As it turns out in the concrete example, we are asking for probabilities from the noncentral t dist in a region of t-values plusminus 2.13 and noncentrality parameters 0.50 and -0.80. These are not at all extreme values. I think it sounds reasonable that "prob1" and "prob2" (using HS's terminology) are around 0.89 and 0.17 [try and plot the densities as function of t, eyeball the relative amount of curve left of the indicated values].
For this reason I am becoming a sceptic now. Either I am not able to make a better implementation than ASA243, or ...please forgive me... the equations given by Julious (and others) are not exact in certain situations (for example: I only see the issue when t2 < nc2), or something entirely else is in play.
Can anyone shed light on this?
Let me add: If the problem lies in ASA243 then I'd love to write one that gets it right, of course.)
Many thanks for any input and best regards,
EM.
Complete thread:
- non-central t, power, R:ElMaestro 2009-10-03 19:05 [🇷 for BE/BA]
- non-central t, power, R: Helmut 2009-10-03 21:45
- non-central t, power, R: ElMaestro 2009-10-04 21:51
- non-central t, power, R: Helmut 2009-10-04 23:58
- non-central t, power, R: ElMaestro 2009-10-04 21:51
- non-central t is not the question d_labes 2009-10-05 12:58
- SAS -> R -> SAS Helmut 2009-10-05 14:53
- Lost in translation d_labes 2009-10-06 14:31
- Lost in translation Helmut 2009-10-06 14:42
- Macacus cynocephalus d_labes 2009-10-06 15:14
- Lost in translation Helmut 2009-10-06 21:46
- Macacus cynocephalus d_labes 2009-10-06 15:14
- Lost in translation Helmut 2009-10-06 14:42
- Lost in translation d_labes 2009-10-06 14:31
- non-central t is not the question ElMaestro 2009-10-05 17:32
- Power of TOST in R d_labes 2009-10-06 16:45
- Power of TOST (R package MBESS) Helmut 2009-10-07 02:18
- Power of TOST (R package MBESS) d_labes 2009-10-07 08:46
- Power of TOST (R package MBESS) Helmut 2009-10-07 12:38
- Power of TOST (R package MBESS) d_labes 2009-10-07 08:46
- The Power at limits d_labes 2009-10-07 11:32
- Power at limits 2 d_labes 2009-10-07 15:15
- package MBESS Helmut 2009-10-07 22:36
- Package MBESS and Power curiosity d_labes 2009-10-15 14:57
- Package MBESS and Power curiosity ElMaestro 2009-10-17 12:26
- Apfelstrudel and EFG 2.01 Helmut 2009-10-17 13:44
- EFG and another power curiosity d_labes 2009-10-20 11:35
- EFG and another power curiosity ElMaestro 2009-11-23 19:03
- EFG and another power curiosity Helmut 2009-11-23 19:17
- EFG and another power curiosity ElMaestro 2009-11-23 19:03
- Package MBESS and Power curiosity ElMaestro 2009-10-17 12:26
- Package MBESS and Power curiosity d_labes 2009-10-15 14:57
- package MBESS Helmut 2009-10-07 22:36
- Power at limits 2 d_labes 2009-10-07 15:15
- Power of TOST (R package MBESS) Helmut 2009-10-07 02:18
- Power of TOST in R d_labes 2009-10-06 16:45
- SAS -> R -> SAS Helmut 2009-10-05 14:53
- non-central t, power, R: yjlee168 2009-10-07 00:08
- R code for power ElMaestro 2009-10-09 11:32
- Belgian beer ? Ohlbe 2009-10-09 11:50
- Another one - R code for power of TROST d_labes 2009-10-09 13:50
- non-central t, power, R: Helmut 2009-10-03 21:45