# Expected P-Value results # R code implementing methods described in Ionita-Laza and Ottman (2011). # The code is useful in comparing the performance of various affected-relative pair designs as well as designs based on unrelated individuals # in the context of association studies with rare variants. # Iuliana Ionita-Laza, June 2011 # Report bugs/comments to ii2135@columbia.edu library(VGAM) l=2 ######################################################################################################################################################### # Function 'EPV_pairs' returns the expected P-value for a study with nnn pairs of affected relatives with exact relationship specified by fi (fi=1/4 for sibs, # fi=1/16 for first cousins, fi=1/64 for second cousins, fi=0 for unrelateds) at a variant with frequency f and genotype relative # risk GRR, for a disease with population prevalence k, and sibling recurrence risk ratio lambda_s. ######################################################################################################################################################### EPV_pairs=function(f,r,k,nnn,fi,lambda_s) { alpha=f*(r-1)/(1+f*(r-1)) #calculate k1 p1a=k/((1-f)*(1-f)*1/r+2*f*(1-f)) k1=2*f*(1-f)*p1a; k2=k-k1; k0=k1/((1-f)*(1-f)+r*2*f*(1-f)); lambda1s=r*k0/k1*2*f*(1-f)*(k0/k1+r*k0/k1)*(2-3*f+f*f)/4+ k0*(1-f)*(1-f)/k1*(k0/k1*(1-f)+r*k0/k1*f); lambda2s=((lambda_s-1)-(k1*k1)/(k*k)*(lambda1s-1))*(k*k)/(k2*k2)+1; if (fi==1/4) { d1=1; lambda=lambda_s; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2s; a1=k1+(1-f)*(1-f)*k2*(r-1)/r; b1=2*f*(1-f)+(1-f)^2/r; x2=a1/b1; x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); nexp=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi)); n=nexp*nnn; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha; n1=n1*nnn; sum_toni=0; for (toni in 0:50) sum_toni=sum_toni+dskellam(toni,n,n1); return (-log10(sum_toni)); } if (fi==1/16) { #first cousins d1=0; lambda=(lambda_s-1)/4+1; lambda1r=(lambda1s-1)/4+1; lambda2r=(lambda2s-1)/4+1; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r; b1=2*f*(1-f)+(1-f)^2/r; x2=a1/b1; x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); nexp=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi)); n=nexp*nnn; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha; n1=n1*nnn; sum_toni=0; for (toni in 0:50) sum_toni=sum_toni+dskellam(toni,n,n1); return(-log10(sum_toni)); } if (fi==1/64) { #second cousins d1=0 lambda=(lambda_s-1)/16+1; lambda1r=(lambda1s-1)/16+1; lambda2r=(lambda2s-1)/16+1; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r; b1=2*f*(1-f)+(1-f)^2/r; x2=a1/b1; x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); nexp=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi)); n=nexp*nnn; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha; n1=n1*nnn; sum_toni=0; for (toni in 0:50) sum_toni=sum_toni+dskellam(toni,n,n1); return(-log10(sum_toni)); } if (fi==0) { #unrelateds d1=0; keff=2; lambda=1; fi=0; lambda1r=1; lambda2r=1; kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); nexp=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi)); n=nexp*nnn; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha; n1=n1*nnn; sum_toni=0; for (toni in 0:50) sum_toni=sum_toni+dskellam(toni,n,n1); return(-log10(sum_toni)); } } ######################################################################################################################################################### # Function 'EPV_one_pairs' returns expected P-value for a study with nnn individuals that have an affected relative with exact relationship specified by # fi (fi=1/4 for sibs, fi=1/16 for first cousins, fi=1/64 for second cousins, fi=0 for unrelateds) at a variant with frequency f and genotype relative # risk GRR, for a disease with prevalence k, and sib recurrence risk ratio lambda_s. ######################################################################################################################################################### EPV_one_pairs=function(f,r,k,nnn,fi,lambda_s) { alpha=f*(r-1)/(1+f*(r-1)) #calculate k1 p1a=k/((1-f)*(1-f)*1/r+2*f*(1-f)) k1=2*f*(1-f)*p1a; k2=k-k1; k0=k1/((1-f)*(1-f)+r*2*f*(1-f)); lambda1s=r*k0/k1*2*f*(1-f)*(k0/k1+r*k0/k1)*(2-3*f+f*f)/4+ k0*(1-f)*(1-f)/k1*(k0/k1*(1-f)+r*k0/k1*f); lambda2s=((lambda_s-1)-(k1*k1)/(k*k)*(lambda1s-1))*(k*k)/(k2*k2)+1; if (fi==1/4) { d1=1; lambda=lambda_s; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2s; a1=k1+(1-f)*(1-f)*k2*(r-1)/r; b1=2*f*(1-f)+(1-f)^2/r; x2=a1/b1; x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); #unrelateds+ n=2*nnn*f; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha*0.5; n1=n1*nnn; sum_toni=0; for (toni in 0:50) sum_toni=sum_toni+dskellam(toni,n,n1); return (-log10(sum_toni)); } if (fi==1/16) { #first cousins d1=0; lambda=(lambda_s-1)/4+1; lambda1r=(lambda1s-1)/4+1; lambda2r=(lambda2s-1)/4+1; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); #unrelateds+ n=2*nnn*f; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha*0.5; n1=n1*nnn; sum_toni=0; for (toni in 0:50) sum_toni=sum_toni+dskellam(toni,n,n1); return (-log10(sum_toni)); } if (fi==1/64) { #second cousins d1=0 lambda=(lambda_s-1)/16+1; lambda1r=(lambda1s-1)/16+1; lambda2r=(lambda2s-1)/16+1; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); #unrelateds+ n=2*nnn*f; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha*0.5; n1=n1*nnn; sum_toni=0; for (toni in 0:50) sum_toni=sum_toni+dskellam(toni,n,n1); return (-log10(sum_toni)); } if (fi==0) { #unrelateds d1=0; keff=2; lambda=1; fi=0; lambda1r=1; lambda2r=1; kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); #unrelateds+ n=2*nnn*f; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha*0.5; n1=n1*nnn; sum_toni=0; for (toni in 0:50) sum_toni=sum_toni+dskellam(toni,n,n1); return (-log10(sum_toni)); } } ######################################################################################################################################################### ######################################################################################################################################################### # The following two functions are needed later on Number_pairs=function(f,r,k,nnn,fi,lambda_s) { alpha=f*(r-1)/(1+f*(r-1)) #calculate k1 p1a=k/((1-f)*(1-f)*1/r+2*f*(1-f)) k1=2*f*(1-f)*p1a; k2=k-k1; k0=k1/((1-f)*(1-f)+r*2*f*(1-f)); lambda1s=r*k0/k1*2*f*(1-f)*(k0/k1+r*k0/k1)*(2-3*f+f*f)/4+ k0*(1-f)*(1-f)/k1*(k0/k1*(1-f)+r*k0/k1*f); lambda2s=((lambda_s-1)-(k1*k1)/(k*k)*(lambda1s-1))*(k*k)/(k2*k2)+1; if (fi==1/4) { d1=1; lambda=lambda_s; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2s; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); nexp=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi)); n=nexp*nnn; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha; n1=n1*nnn; vn=c(n,n1); return (vn); } if (fi==1/16) { #first cousins d1=0; lambda=(lambda_s-1)/4+1; lambda1r=(lambda1s-1)/4+1; lambda2r=(lambda2s-1)/4+1; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); nexp=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi)); n=nexp*nnn; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha; n1=n1*nnn; vn=c(n,n1); return (vn); } if (fi==1/64) { #second cousins d1=0 lambda=(lambda_s-1)/16+1; lambda1r=(lambda1s-1)/16+1; lambda2r=(lambda2s-1)/16+1; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); nexp=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi)); n=nexp*nnn; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha; n1=n1*nnn; vn=c(n,n1); return (vn); } if (fi==0) { #unrelateds d1=0; keff=2; lambda=1; fi=0; lambda1r=1; lambda2r=1; kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); nexp=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi)); n=nexp*nnn; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*keff*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha; n1=n1*nnn; vn=c(n,n1); return (vn); } } ######################################################################################################################################################### ######################################################################################################################################################### Number_one_pairs=function(f,r,k,nnn,fi,lambda_s) { alpha=f*(r-1)/(1+f*(r-1)) #calculate k1 p1a=k/((1-f)*(1-f)*1/r+2*f*(1-f)) k1=2*f*(1-f)*p1a; k2=k-k1; k0=k1/((1-f)*(1-f)+r*2*f*(1-f)); lambda1s=r*k0/k1*2*f*(1-f)*(k0/k1+r*k0/k1)*(2-3*f+f*f)/4+ k0*(1-f)*(1-f)/k1*(k0/k1*(1-f)+r*k0/k1*f); lambda2s=((lambda_s-1)-(k1*k1)/(k*k)*(lambda1s-1))*(k*k)/(k2*k2)+1; if (fi==1/4) { d1=1; lambda=lambda_s; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2s; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); #unrelateds+ n=2*nnn*f; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha*0.5; n1=n1*nnn; vn=c(n,n1); return (vn); } if (fi==1/16) { #first cousins d1=0; lambda=(lambda_s-1)/4+1; lambda1r=(lambda1s-1)/4+1; lambda2r=(lambda2s-1)/4+1; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); #unrelateds+ n=2*nnn*f; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha*0.5; n1=n1*nnn; vn=c(n,n1); return (vn); } if (fi==1/64) { #second cousins d1=0 lambda=(lambda_s-1)/16+1; lambda1r=(lambda1s-1)/16+1; lambda2r=(lambda2s-1)/16+1; keff=log(4*f*fi+4*f*f*(1-4*fi+4*d1*fi*fi))/log(2*f); kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); #unrelateds+ n=2*nnn*f; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha*0.5; n1=n1*nnn; vn=c(n,n1); return (vn); } if (fi==0) { #unrelateds d1=0; keff=2; lambda=1; fi=0; lambda1r=1; lambda2r=1; kr=k*lambda; k2r=k2*lambda2r; a1=k1+(1-f)*(1-f)*k2*(r-1)/r b1=2*f*(1-f)+(1-f)^2/r x2=a1/b1 x1=(k1-2*f*(1-f)*x2)/(1-f)^2; beta=(x2*x2+2*x2*(1-x2)*k2+(1-x2)^2*k2*k2r)/(k*kr); alpha=(1-(1-x1)*(1-k2)-(1-x2)*(1-k2)+(1-x1)*(1-x2)*(1-2*k2+k2*k2r))/(k*kr); #unrelateds+ n=2*nnn*f; n1=4*f*(fi+f*(1-4*fi+4*d1*fi*fi))*beta+4*f*(1-f)*(1-f)*(1-2*fi-f*(1-4*fi+4*d1*fi*fi))*alpha*0.5; n1=n1*nnn; vn=c(n,n1); return (vn); } } ######################################################################################################################################################### # Function 'EPV_pairs_aggregate' returns expected P-value for a study with nnn pairs of affected relatives with exact relationship specified by fi # (fi=1/4 for sibs, fi=1/16 for first cousins, fi=1/64 for second cousins, fi=0 for unrelateds) at a locus with nd variants with frequencies # given in vector f and genotype relative risk in vector r for a disease with prevalence k, and sib recurrence risk ratio lambda_s ######################################################################################################################################################### EPV_pairs_aggregate=function(f,r,nd,k,nnn,fi,lambda_s) { n_t=0; n1_t=0; for (i in 1:nd) { n_t=n_t+Number_pairs(f[i],r[i],k,nnn,fi,lambda_s)[1]; n1_t=n1_t+Number_pairs(f[i],r[i],k,nnn,fi,lambda_s)[2]; } sum_toni=0; for (toni in 0:50) sum_toni=sum_toni+dskellam(toni,n_t,n1_t); return (-log10(sum_toni)); } ######################################################################################################################################################### # Function 'EPV_one_pairs_aggregate' returns expected P-value for a study with nnn affected individuals that have an affected relative with exact # relationship specified by fi (fi=1/4 for sibs, fi=1/16 for first cousins, fi=1/64 for second cousins, fi=0 for unrelateds) at a locus with nd # variants with frequencies given in vector f and genotype relative in vector r for a disease with prevalence k, and sib recurrence risk ratio lambda_s ######################################################################################################################################################### EPV_one_pairs_aggregate=function(f,r,nd,k,nnn,fi,lambda_s) { n_t=0; n1_t=0; for (i in 1:nd) { n_t=n_t+Number_one_pairs(f[i],r[i],k,nnn,fi,lambda_s)[1]; n1_t=n1_t+Number_one_pairs(f[i],r[i],k,nnn,fi,lambda_s)[2]; } sum_toni=0; for (toni in 0:50) sum_toni=sum_toni+dskellam(toni,n_t,n1_t); return (-log10(sum_toni)); } ######################################################################################################################################################### ######################################################################################################################################################### ######################################################################################################################################################### ######################################################################################################################################################### # Examples # To calculate the expected P-value at a variant with frequency f=0.01 and GRR=1.2 for a study with 1000 affected sib-pairs, # for a disease with population prevalence k=0.03 and sib recurrence risk ratio lambda_s=2 type: # > EPV_pairs(0.01,1.2,0.03,1000,1/4,2) # Similarly, for a study with 2000 affected individuals with an affected sibling, type: # > EPV_one_pairs(0.01,1.2,0.03,2000,1/4,2) # For a locus with 3 disease variants and 3 random variants with f=(0.001,0.005,0.01,0.001,0.001,0.007) and GRR=(2,2,2,1,1,1) and # for a study with 1000 affected sib-pairs, for a disease with population prevalence k=0.03 and sib recurrence risk ratio lambda_s=2 type: # > f=c(0.001,0.005,0.01,0.001,0.001,0.007) # > GRR=c(2,2,2,1,1,1) # > EPV_pairs_aggregate(f,GRR,6,0.03,1000,1/4,2)