Multiple comparison adjustment where alpha changes across tests

[sadf...@gmail.com]

Hello:

prob(at least one false pos) = 1 - (1 - alpha)^k

where

k = number of tests

How would this formula be adjusted if alpha changed depending on the test?

Example:

5 tests are performed:
4 of the 5 tests are set at alpha=.05
1 of the 5 tests is set at alpha=.01

[Ray Koopman]

(1 - alpha_1)*(1 - alpha_2)* ... *(1 - alpha_k)

[sadf...@gmail.com]

thanks. I read that this formula is derived from the binomial pdf. How is it connected to the binomial pdf?

also, how could we adjust the formula above to take into account non-independence among tests?--what if we knew the correlation between the tests (e.g., rho(test1,test2) = .10; rho(test1,test3)=.15, etc.)?

thanks for your help.

[Ray Koopman]

On Jan 10, 2:44 am, sadf...@gmail.com wrote:
> On Thursday, January 10, 2013 12:56:19 AM UTC-5, Ray Koopman wrote:

>> On Jan 9, 3:30 pm, sadf...@gmail.com wrote:
>>> Hello:
>>> prob(at least one false pos) = 1 - (1 - alpha)^k
>>> where
>>> k = number of tests
>>> How would this formula be adjusted if alpha changed depending on
>>> the test?
>>> Example:
>>> 5 tests are performed:
>>> 4 of the 5 tests are set at alpha=.05
>>> 1 of the 5 tests is set at alpha=.01
>>
>> prob(at least one false pos) = 1 -
>> (1 - alpha_1)*(1 - alpha_2)* ... *(1 - alpha_k)
>
> thanks. I read that this formula is derived from the binomial pdf.
> How is it connected to the binomial pdf?

Look up the binomial distribution.
If X is Binomial(k,p) then Prob(X = 0) = (1-p)^k,
from which it follows that Prob(X > 0) = 1 - (1-p)^k.

>
> also, how could we adjust the formula above to take into account non
> -independence among tests?--what if we knew the correlation between
> the tests (e.g., rho(test1,test2) = .10; rho(test1,test3)=.15, etc.)?

There is no exact general answer. The best simple answer is what is
usually referred to as the Bonferroni upper bound:

prob(at least one false pos) <= alpha_1 + alpha_2 + ... + alpha_k.

There have been many improvements suggested, but all of them are more
complex. Look up multiple comparisons.

[sadf...@gmail.com]

very helpful. thank you. i have another question. what if we wanted to know the probability of the following:

1. exactly two false positives for k independent tests
2. at least two false positives for k independent tests

assuming that the alphas stay the same across the tests and then assuming that the alphas change.

i've examined the binomial distribution, but I'm still struggling. this is not homework. just trying to learn.

[sadf...@gmail.com]

actually, i think i just figured out my question WHEN alpha stays the same. if the alphas are the same then i can apply the formula:

[n!/(n-k)!]*(p^k)*[(1-p)^(n-k)]

where
n is the total number of tests
k is the number of false positive tests
p is alpha

So, if we are solving for prob(k=0), then [n!/(n-k)!]*(p^k) = 1, reducing to:

(1 - p)^n

then it makes sense that 1 - (1 - p)^n = prob(k>0)

if we are solving for prob(k=1), prob(k=2) etc., then [n!/(n-k)!]*(p^k) does not cancel automatically. then i just have to incorporate the factorial term.

but not sure how to solve, for example, prob(k=2), when p varies???

[Ray Koopman]

On Jan 10, 7:27 pm, sadfjl...@gmail.com wrote:
> On Thursday, January 10, 2013 9:44:05 PM UTC-5, sadf...@gmail.com wrote:
>> On Thursday, January 10, 2013 8:29:28 PM UTC-5, Ray Koopman wrote:

>>> On Jan 10, 2:44 am, sadfjl...@gmail.com wrote:
>>>> On Thursday, January 10, 2013 12:56:19 AM UTC-5, Ray Koopman wrote:

You dropped the k! in the denominator of the binomial coefficient,
n_Choose_k = n! /( k! * (n-k)! ).

>
> where
> n is the total number of tests
> k is the number of false positive tests
> p is alpha
>
> So, if we are solving for prob(k=0), then [n!/(n-k)!]*(p^k) = 1,
> reducing to:
>
> (1 - p)^n
>
> then it makes sense that 1 - (1 - p)^n = prob(k>0)
>
> if we are solving for prob(k=1), prob(k=2) etc.,
> then [n!/(n-k)!]*(p^k) does not cancel automatically.
> then i just have to incorporate the factorial term.
>
> but not sure how to solve, for example, prob(k=2), when p varies???

When p varies, you're no longer dealing with the ordinary
binomial distribution. Then it's called the compound binomial,
and the computations are messy. You have to sum the separate
probabilities of the n_Choose_k possible patterns with k false
positives. Here's an example for n = 5, k = 2.

( p1 ) * ( p2 ) * (1-p3) * (1-p4) * (1-p5) +
( p1 ) * (1-p2) * ( p3 ) * (1-p4) * (1-p5) +
( p1 ) * (1-p2) * (1-p3) * ( p4 ) * (1-p5) +
( p1 ) * (1-p2) * (1-p3) * (1-p4) * ( p5 ) +
(1-p1) * ( p2 ) * ( p3 ) * (1-p4) * (1-p5) +
(1-p1) * ( p2 ) * (1-p3) * ( p4 ) * (1-p5) +
(1-p1) * ( p2 ) * (1-p3) * (1-p4) * ( p5 ) +
(1-p1) * (1-p2) * ( p3 ) * ( p4 ) * (1-p5) +
(1-p1) * (1-p2) * ( p3 ) * (1-p4) * ( p5 ) +
(1-p1) * (1-p2) * (1-p3) * ( p4 ) * ( p5 )

If all the pj were equal then all the terms in the above sum would be
equal and it would simplify to the usual 5_Choose_2 * p^2 * (1-p)^3.

[sadf...@gmail.com]

thank you very much. this clears everything up for me. Just one last question re: the solution you provided in your initial response (which makes sense, of course):

prob(at least one false pos) = 1 - ( 1 - alpha_1 ) * ( 1 - alpha_2) * ... * (1 - alpha_k)

is that from the standard binomial (and simply change the equation after canceling the out the factorial) or did you arrive at this equation from the compound binomial?

thanks again!

[sadf...@gmail.com]

actually another couple questions if its ok...

what if we wanted to solve

prob(k=1 or k=2)=? for both scenarios (constant p and varying p).

Also,

prob(k=1 or k= 2) for constant p but varying n (# of trials). in other words, what's the

prob(k=1 for n=3 OR k=2 for n=4) =?

thanks again and sorry for all the questions! This is really interesting tho! :-/

[Ray Koopman]

compound binomial

[Ray Koopman]

On Jan 11, 5:55 am, sadfjl...@gmail.com wrote:
> what if we wanted to solve
>
> prob(k=1 or k=2)=? for both scenarios (constant p and varying p).

prob(k=1 or k=2) = prob(k=1) + prob(k=2).

>
> Also,
>
> prob(k=1 or k= 2) for constant p but varying n (# of trials).
> in other words, what's the
>
> prob(k=1 for n=3 OR k=2 for n=4) =?

There is no solution unless you say how the choice of n is made.

[sadf...@gmail.com]

Hi,

thank you for clearing up so many of my questions. i can see how I'll be able to use this distribution quite often. problem is that i can't find the pdf and log likelihood of this distribution anywhere (textbooks, internet). Would you mind sharing them with me?

thanks so much for your time.

[Ray Koopman]

X is a Compound Binomial variable with parameters p_1,...,p_n
if it is the sum of n independent Bernoulli variables whose
parameters are p_1,...,p_n. (If p_1 = ... = p_n = p then
X is an ordinary Binomial variable with parameters n,p.)
Prob(X = k) = the coefficient of s^k in the expansion of
product{j = 1,...,n: (1-p_j + s*p_j) }.

[sadf...@gmail.com]

let's see if i got this straight with a very simple example...

Ordinary Binomial pdf:

Prob(X = k) = (n!/(k!*(n-k)!))*(p^k)*(1-p)^(n-k)

where
n is the total number of tests
k is the number of false positive tests
p is alpha

Example:

n = 2
k = 1
p = .05

Prob = (2!/(1!*(2-1)!))*(.05^1)*(1-.05)^(2-1) = .095

Compound Binomial pdf (based on your description):

Prob(X = k) = the coefficient of s^k in the expansion of the product{j = 1,...,n: (1-p_j + s*p_j)}

Example:

n = 2
k = 1
p_1 = .05
p_2 = .10

Prob = (.05)*(1-.10) + (.10)*(1-.05) = .14

this doesn't seem right to me. can you tell me where I've gone off the rail? By the way, what does "s" stand for in the compound binomial pdf? i think that might be tripping me up. sorry--i think i got mixed up with the notation.

thanks!

[Ray Koopman]

Let qj = 1 - pj. Then expand (q1 + s*p1)*(q2 + s*p2)

= q1*q2 + (q1*p2 + p1*q2)*s + p1*p2*s^2.

s is a dummy variable whose function is to group the products
of the p's and q's on the basis of how many p's they contain,
which is both the exponent of s and the value of x.

Pr(x = 0) = q1*q2
Pr(x = 1) = q1*p2 + p1*q2
Pr(x = 2) = p1*p2

If s = 1 then the numeric value of the whole expression is 1,
which is another way of saying that the probabilities sum to 1.

[sadf...@gmail.com]

thanks! now i'm starting to make the connection! not to prolong this discussion, but I'm curious if you know how one would go about taking the log of the likelihood function of the compound binomial distribution in order to employ MLE.

For the standard binomial distribution, the pdf is:

p(y,mu) = c(n,y) * mu^y * (1-mu)^(n-y)

we can maximize the log of the pdf above or, as many stats software packages, we can omit the coefficient (constant) and maximize the log of:

f(p(y,mu)) = mu^y * (1-mu)^(n-y)

Using the NLMIXED procedure in SAS, we could maximize the log of p(y,mu) like this:

proc nlmixed data={dataset};
eta = b0;
p = exp(eta)/(1 + exp(eta));
prob = comb(n,y) * (p**y) * (1-p)**(n-y);
loglike = log(prob);
model y ~ general(loglike);
run;

In the example above, the linear predictor (eta) just includes an intercept (b0), but if we wanted, we could add other predictors. Regardless, we could also maximize f(p(y,mu)), like:

proc nlmixed data={dataset};
eta = b0;
p = exp(eta)/(1 + exp(eta));
prob = (p**y) * (1-p)**(n-y);
loglike = log(prob);
model y ~ general(loglike);
run;

both will yield the same estimate of mu since the coefficient is simply a constant. anyway, point is that i'm not sure how to obtain the log the pdf of the compound binomial distribution.

I'm not sure if you use SAS or some other program, but would you mind sharing with me equations/(code for whatever program you use) in order to take the log of the pdf of the compound binomial distribution from which one would employ MLE. Any hints/tips would be great!

i understand if you don't have time to respond anymore on this topic. i must say you're a great teacher, even over the internet. i hope you do this as a full time job!

[Ray Koopman]

On Jan 17, 9:23 am, sadfjl...@gmail.com wrote:
>
> thanks! now i'm starting to make the connection! not to prolong this discussion, but I'm curious if you know how one would go about taking the log of the likelihood function of the compound binomial distribution in order to employ MLE.
>
> For the standard binomial distribution, the pdf is:
>
> p(y,mu) = c(n,y) * mu^y * (1-mu)^(n-y)
>
> we can maximize the log of the pdf above or, as many stats software packages, we can omit the coefficient (constant) and maximize the log of:
>
> f(p(y,mu)) = mu^y * (1-mu)^(n-y)
>
> Using the NLMIXED procedure in SAS, we could maximize the log of p(y,mu) like this:
>
> proc nlmixed data={dataset};
> eta = b0;
> p = exp(eta)/(1 + exp(eta));
> prob = comb(n,y) * (p**y) * (1-p)**(n-y);
> loglike = log(prob);
> model y ~ general(loglike);
> run;
>
> In the example above, the linear predictor (eta) just includes an intercept (b0), but if we wanted, we could add other predictors. Regardless, we could also maximize f(p(y,mu)), like:
>
> proc nlmixed data={dataset};
> eta = b0;
> p = exp(eta)/(1 + exp(eta));
> prob = (p**y) * (1-p)**(n-y);
> loglike = log(prob);
> model y ~ general(loglike);
> run;
>
> both will yield the same estimate of mu since the coefficient is simply a constant. anyway, point is that i'm not sure how to obtain the log the pdf of the compound binomial distribution.
>
> I'm not sure if you use SAS or some other program, but would you mind sharing with me equations/(code for whatever program you use) in order to take the log of the pdf of the compound binomial distribution from which one would employ MLE. Any hints/tips would be great!

You need to expand the probabilities -- by hand, or (preferably) by a
computer algebra system. I use Mathematica. Here is sample code and
the output it gives for n = 3:

CoefficientList[Product[q[j]+s*p[j],{j,n}],s] /. q[j_] -> 1-p[j]

{(1-p[1])*(1-p[2])*(1-p[3]),
p[1]*(1-p[2])*(1-p[3]) + (1-p[1])*p[2]*(1-p[3]) +
(1-p[1])*(1-p[2])*p[3],
p[1]*p[2]*(1-p[3]) + p[1]*(1-p[2])*p[3] + (1-p[1])*p[2] p[3],
p[1]*p[2]*p[3]}

Call those expressions pr[0],pr[1],pr[2],pr[3], because they are the
probabilities of y = 0,1,2,3, expressed as functions of the 3 unknowns
p[1],p[2],p[3]. Then loglike = log(pr[y]). You would also need
something like p[j] = 1/(1 + e^(-eta[j])), j = 1,2,3, possibly with
each eta[j] being a different linear function of some independent
variables. I don't know SAS, so I can't help you with the details of
communicating the model to nlmixed. Maybe someone else on the list
will help out here.

[sadf...@gmail.com]

thank you! I'm not familiar enough with nlmixed to not have to expand the probabilities myself. with that said, now that you've shown me what needs to be done, I think I know how to write the nlmixed code if I expand the probabilities myself. but, first, i want to make sure I understand what the dataset would need to look like:

id test1 test2 test3 y
1 0 0 0 0
2 1 1 1 3
3 0 1 0 1
4 1 0 1 2
5 0 0 0 0
.
.
.

where for each test

0 = not a false positive
1 = false positive

and for y

0 = no false positives
1 = exactly 1 false positive
2 = exactly 2 false positives
3 = exactly 3 false positives

Does the above look right?

anyway, the nlmixed code would look like this, i think:

proc nlmixed data = {dataset};

eta1 = beta1_0;
eta2 = beta2_0;
eta3 = beta3_0;

p1 = 1 / (1 + exp(-eta1));
p2 = 1 / (1 + exp(-eta2));
p3 = 1 / (1 + exp(-eta3));

pr0 = (1-p1)*(1-p2)*(1-p3);
pr1 = p1*(1-p2)*(1-p3) + (1-p1)*p2*(1-p3) + (1-p1)*(1-p2)*p3;
pr2 = p1*p2*(1-p3) + p1*(1-p2)*p3 + (1-p1)*p2*p3;
pr3 = p1*p2*p3;

if y = 0 then loglike = log(pr0); else
if y = 1 then loglike = log(pr1); else
if y = 2 then loglike = log(pr2); else
if y = 3 then loglike = log(pr3);

model y ~ general(loglike);

run;

yea, i think above is correct. PLEASE let me know if you disagree. otherwise, thanks for all your help. this has been really fun and a great learning experience. Cool. :-)

[Ray Koopman]

I don't know how nlmixed knows which variables' values you want it to
choose so as to maximize the sum of the loglike's, but otherwise it
looks OK. However, for a "location only" model such as you are using,
the MLE of each pj is simply the corresponding proportion of 1's in
testj.