A Bayesian decision procedure for testing multiple hypotheses in DNA microarray experiments

Appendix A

Proof of Proposition 1

First we show that for fixed values of the costs for false negatives, C_0i, i=1, …, N, the posterior expected loss for the Bayes rule that rejects k hypotheses, is a function of the cutoff p_k or equivalently a function of k, the number of rejected null hypotheses.

In fact, setting p_k=C_0i/(C_0i+C_1i), we can write as follows,

where is the posterior Bayes action by which the kth null hypothesis will be rejected. Then, for fixed values of C_0i, i=1, …, N, the posterior expected loss for the Bayes rule, is a function of the cutoff p_k or equivalently a function of k, the number of rejected null hypotheses.

Finally we show that is a decreasing function on k.

Let k₁ and k₂ be such that k₁<k₂; then and and

Appendix B

In fact, if C_0i=C for i=1, …, N, then the second term of (19) can be written as

where p=Pr(H_0i=0|p) is the prior probability that each null hypothesis is true. Then, substituting (20) in (19), we immediately obtain (13).

Appendix C

Matlab code for simulated data

 1. Codes of Matlab for function M–file for function to maximize:

1 function EB = EB (c,t,al,bet,a,b,m,N)

2 p= ;phi= ;mu=zeros(1,N);EB=0;burnin= ; iters= ;

3 for iter=1:burnin+iters

4  prob0=p*exp(-phi*(t.^2)/2);

5  prob1=prob0+(1-p)*exp(-phi*((t-mu).^2)/2);

6  uni=rand(1,N).*prob1;z=zeros(1,N);z=uni>prob0;

7  N1=sum(z);N0=N-N1;I1=find(z);

8  p=betarnd(al+N0,bet+N1);

9  bast=b+sum(t(z==0).^2)+sum((t(I1)-mu(I1)).^2)+sum(c*(mu-m).^2);

10  phi=gamrnd((a+2*N)/2,2/bast);

11  mu(I1)=normrnd((c*m(I1)+t(I1))./(c+1),1./sqrt((c+1)*phi));

23  mu(z==0)=normrnd(m(z==0),1./sqrt(c*phi));

13  if iter>burnin;EB=EB+prod((p*(phi^0.5)*exp(-phi*(t.^2)/2))+((1-14 p)*(phi^0.5)* exp(-phi*((t-mu).^2)/2)));end

15 end

16 EB=-EB/iters

 2. Create an M–file to obtain the maximum of the function EB using the code “fminbnd”

 3. Given the previously estimated value of c, next sentences estimate the rest of the parameters and measures used in the approach.

17 clear;seed=;rand(’state’,seed);randn(’state’,seed);

18 N=;n=;sig=;ptrue=;phitrue=n*(sig)^-2;mutrue=linspace( , ,N);

19 al=;bet=;a=;b=;m=zeros(1,N);c= *ones(1,N);

20 meds=mutrue;uni=rand(1,N);ztrue=uni>ptrue;meds(uni<ptrue)=0;

21 for i=1:n;x(i,:)=normrnd(meds,sig);end

22 t=mean(x);s=std(x);

23 p=; phi=; mu=zeros(1,N);

24 Ez=zeros(1,N);Emu1=zeros(1,N);Ep=0;Ephi=0;

25 fid = fopen(’pphi.txt’,’wt’);fid2 = fopen(’tzmu.txt’,’wt’);

26 fid3 = fopen(’mpphi.txt’,’wt’);burnin=;iters=;

27 for iter=1:burnin+iters

28  prob0=p*exp(-phi*(t.^2)/2);

29  prob1=prob0+(1-p)*exp(-phi*((t-mu).^2)/2);

30  uni=rand(1,N).*prob1;z=zeros(1,N);z=uni>prob0;

31  if iter>burnin;Ez=Ez+z;end

32  N1=sum(z);N0=N-N1;I1=find(z);p=betarnd(al+N0,bet+N1);

33  if iter>burnin;Ep=Ep+p;mp=Ep/(iter-burnin);end

34  bast=b+sum(t(z==0).^2)+sum((t(I1)-mu(I1)).^2)+sum(c.*(mu-m).^2);

35  phi=gamrnd((a+2*N)/2,2/bast);

36  if iter>burnin;Ephi=Ephi+phi;mphi=Ephi/(iter-burnin);end

37  if iter>burnin;fprintf(fid3,’%6.4f %6.4f/n’,[mp; mphi]);end

38  mu(I1)=normrnd((c(I1).*m(I1)+t(I1))./(c(I1)+1),1./sqrt((c(I1)+1)*phi));

39  mu(z==0)=normrnd(m(z==0),1./sqrt(c(z==0)*phi));

40  if iter>burnin;Emu1(I1)=Emu1(I1)+mu(I1);end

41  if iter>burnin;fprintf(fid,’%6.4f %6.4f/n’,[p; phi]);end

42 end

43 Emu1=Emu1./Ez;Ez=Ez/iters;load pphi.txt;load mpphi.txt

44 pest=[ 1-sum(ztrue)/N ptrue mean(pphi(:,1))]

45 phiest=[phitrue mean(pphi(:,2))]

46 zest=[ztrue;Ez]’;muest=[meds;Emu1]’;tzmu=[t; ztrue;Ez;meds;Emu1]’;

47 fprintf(fid2,’%6.4f %6.4f %6.4f %6.4f %6.4f/n’,[t;ztrue;Ez;meds;Emu1]);

48 fclose(’all’);

49 subplot(2,2,1),plot(pphi(:,1)),title(’p’)

50 subplot(2,2,2),plot(pphi(:,2)),title(’phi’)

51 subplot(2,2,3),plot(mpphi(:,1)),title(’mp’)

52 subplot(2,2,4),plot(mpphi(:,2)),title(’mphi’)

53 RB=(sum(Ez>.5)*100)/N

54 R1=round(N*(1-mean(pphi(:,1))));

55 N1=(R1*100)/N

56 EzO=sort(Ez,’descend’);pN1=1-EzO(R1)

57 V=cumsum(1-EzO);F=cumsum(EzO);T=sum(Ez)-F;FPr=V./V(N);

58 FNr=T./sum(Ez);R=(1:1:N);I=V./R;II=T./(N-R);

59 for i=1:N-1;EI(i)=I(i);EII(i)=II(i);end

60 rFDR=[EI,V(N)/N];rFNR=[EII,0];

61 FPrB=FPr((RB*N)/100)

62 FNrB=FNr((RB*N)/100)

63 rFDRB=rFDR((RB*N)/100)

64 rFNRB=rFNR((RB*N)/100)

65 FPrFNH=FPr(R1)

66 FNrFNH= FNr(R1)

67 rFDRFNH=rFDR(R1)

68 rFNRFNH=rFNR(R1)

• For real data change the following lines

• line 17 by: clear;x = xlsread(’File name.xls’); %This file must contain a column with mean differences

• line 18 by: N= ;

• lines 20, 21 and 22 by: t=x’;

• line 44 by: pest=mean(pphi(:,1))

• line 45 by: phiest=mean(pphi(:,2))

• line 46 by: tzmu=[t;Ez;Emu1]’;

• line 47 by: fprintf(fid2,’%6.4f %6.4f %6.4f/n’,[t;Ez;Emu1]);