Nothing
prop.comb <-
function(x,n,p=NULL,a=NULL,alternative = c("two.sided", "less",
"greater"), conf.level = 0.95, coverage=FALSE, nrep=1000)
{
noa=missing(a)
a0=a
# -------------------------------------------
if (length(x) != length(n))
stop("'x' and 'n' must have the same length")
if (any(n <= 0))
stop("elements of 'n' must be positive")
if (any(x < 0))
stop("elements of 'x' must be nonnegative")
if (any(x > n))
stop("elements of 'x' must not be greater than those of 'n'")
xr <- round(x)
if (any(is.na(x) | (x < 0)) || max(abs(x - xr)) > 1e-07)
stop("'x' must be nonnegative and integer")
x=xr
nr <- round(n)
if (any(is.na(n) | (n < 1)) || max(abs(n - nr)) > 1e-07)
stop("'n' must be a positive integer")
n=nr
if (!missing(a)) {if (length(a) != length(n))
stop("'a' and 'n' must have the same length")}
#---------------------------------------------
########################################
# AUXILIAR FUNCTIONS
########################################
#### One proportion: test1.MA + ic1.MA
#####################
test1.MA=function(x,n,p0=0.5,e=5,alternative = c("two.sided", "less", "greater")) {
c=1.0/(2.0*n) ; pp=x/n
dif=pp-p0; part1= dif*dif*n
part2=p0*(1-p0) ;ze0=part1/part2
difc=abs(dif)
if (difc>c) {ze0c=n*(difc-c)*(difc-c)/part2} else {ze0c=0.0}
b1=sqrt((x+0.5)/(n+1.0)); b2=sqrt(p0)
fact=asin(b1)-asin(b2)
a1=4.0*(n+1.0)*fact*fact
fact1=(x+2.0-(n+4.0)*p0); fact2=(x+2.0)*(n-x+2.0);
zw2=(fact1*fact1)*(n+4.0)/fact2;
cc1=sqrt(ze0); cc2=sqrt(ze0c); cc3=sqrt(a1); cc4=sqrt(zw2);
alternative <- match.arg(alternative)
if(pp<p0){cc1=-cc1;cc2=-cc2;cc3=-cc3;cc4=-cc4} # eran todos positivos
p_ze0=switch(alternative,two.sided = 2.0*pnorm(-abs(cc1)),less = pnorm(cc1),greater = 1-pnorm(cc1))
p_ze0c=switch(alternative,two.sided = 2.0*pnorm(-abs(cc2)),less = pnorm(cc2),greater = 1-pnorm(cc2))
p_a1=switch(alternative,two.sided = 2.0*pnorm(-abs(cc3)),less = pnorm(cc3),greater = 1-pnorm(cc3))
p_zw2=switch(alternative,two.sided = 2.0*pnorm(-abs(cc4)),less = pnorm(cc4),greater = 1-pnorm(cc4))
z=numeric(4)
pvalue=numeric(4)
z[1]=sqrt(ze0); pvalue[1]=p_ze0 # Score (without cc)
z[2]=sqrt(ze0c); pvalue[2]=p_ze0c # Score (with cc)
z[3]=sqrt(a1); pvalue[3]=p_a1 # Adjusted Arc Sine (a)
z[4]=sqrt(zw2);pvalue[4]=p_zw2 # Adjusted Wald
result=cbind(z,pvalue)
rownames(result)=c("Score (without cc)","Score (with cc)",
"Adjusted Arc Sine", "Adjusted Wald")
colnames(result)=c("z value", "p-value")
result}
ic1.MA=function(x,n,e=5) {
c=1.0/(2.0*n)
e1=1.0-(e/100.0)/2.0; Z=qnorm(e1); pp=x/n
part1=n/(n+Z*Z);
part2=(Z*Z)/(4.0*n*n)+pp*(1-pp)/n;
part3=pp+Z*Z/(2.0*n);
l1ze0=part1*(part3-Z*sqrt(part2));
l2ze0=part1*(part3+Z*sqrt(part2));
part2ccn=(Z*Z)/(4.0*n*n)+(pp-c)*(1-pp+c)/n;
part2ccp=(Z*Z)/(4.0*n*n)+(pp+c)*(1-pp-c)/n;
l1ze0c=part1*(part3-c-Z*sqrt(part2ccn));
l2ze0c=part1*(part3+c+Z*sqrt(part2ccp));
b1=sqrt((x+0.5)/(n+1.0)); b2=Z/(2.0*sqrt(n+1.0))
factl1=asin(b1)-b2; factl2=asin(b1)+b2
l1a1=sin(factl1)*sin(factl1); l2a1=sin(factl2)*sin(factl2)
fact1=(x+2.0)*(n+2.0-x)/(n+4.0);
fact2=(x+2.0)-Z*sqrt(fact1);
fact3=(x+2.0)+Z*sqrt(fact1);
l1zw2=fact2/(n+4.0);
l2zw2=fact3/(n+4.0);
part11=(n+Z*Z*Z*Z/53)/(n+Z*Z);
part21=Z/(n+Z*Z);
raiz1=n*pp*(1-pp)+Z*Z/4;
l1y0=0.5+part11*(pp-0.5)-part21*sqrt(raiz1);
l2y0=0.5+part11*(pp-0.5)+part21*sqrt(raiz1);
if (l1ze0<0.0) l1ze0=0.0; if (l2ze0>1.0) l2ze0=1.0;
if (l1ze0c<0.0) l1ze0c=0.0; if (l2ze0c>1.0) l2ze0c=1.0;
if (l1a1<0.0) l1a1=0.0; if (l2a1>1.0) l2a1=1.0;
if (l1zw2<0.0) l1zw2=0.0; if (l2zw2>1.0) l2zw2=1.0;
if (l1y0<0.0) l1y0=0.0; if (l2y0>1.0) l2y0=1.0;
result=matrix(nrow=5,ncol=2)
result[1,]=c(l1ze0,l2ze0)
result[2,]=c(l1ze0c,l2ze0c)
result[3,]=c(l1a1,l2a1)
result[4,]=c(l1zw2,l2zw2)
result[5,]=c(l1y0,l2y0)
rownames(result)=c("Score (without cc)","Score (with cc)",
"Adjusted Arc Sine", "Adjusted Wald", "Modified Score" )
colnames(result)=c("lower limit", "upper limit")
result
}
prop.test.MA=function(x,n,p0=0.5,conf.level = 0.95,
alternative = c("two.sided", "less", "greater")) {
e=100*(1-conf.level)
A=test1.MA(x,n,p0,e,alternative)
alternative <- match.arg(alternative)
B=switch(alternative,two.sided = ic1.MA(x,n,e),less = ic1.MA(x,n,2*e),greater = ic1.MA(x,n,2*e))
B[,1]=switch(alternative,two.sided = B[,1],less = 0,greater=B[,1])
B[,2]=switch(alternative,two.sided = B[,2],less = B[,2],greater=1)
A=rbind(A,c(NA,NA))
result=cbind(B,A)
# C=binom.test(x,n,p=p0,alternative=alternative)
# C=c(as.numeric(C$conf.int),NA,C$p.value)
# result=rbind(result,C)
# rownames(result)[6]="Exact binomial test"
list(estimate=x/n,inference=result,alternative=alternative,p0=p0)}
#############################################
# ESTATISTIC BASED ON ARCSINE TRANSFORMATION
#############################################
As=function(delta,x2,n2,a,fact1,e1) {
p_vero1=Pvero_dif(delta=delta*sign(a[2]),x=x2,n=n2)
p_vero2=p_vero1 +delta/sign(a[2])
p_vero2[p_vero2<0.0]=0; p_vero2[p_vero2>1.0]=1
fact2=asin(sqrt(p_vero2))-asin(sqrt(p_vero1))
qnorm(e1)^2-4.0*n2[1]*n2[2]*(fact1-fact2)*(fact1-fact2)/(n2[1]+n2[2])}
AM=function(delta,x,n2,a,e1) {
p_vero1=Pvero_dif(delta=delta*sign(a[2]),x=x,n=n2)
p_vero2=p_vero1+delta/sign(a[2])
p_vero2[p_vero2<0.0]=0; p_vero2[p_vero2>1.0]=1
b1=sqrt(p_vero1); b2=sqrt(p_vero2)
b3=sqrt(x[1]/n2[1]); b4=sqrt(x[2]/n2[2])
qnorm(e1)^2-4.0*(n2[1]*(asin(b1)-asin(b3))^2+n2[2]*(asin(b2)-asin(b4))^2)}
##############################################
#### Difference of proportions: test2.MA + ic2.MA
###############################
ic2.MA=function(x,n,e=5,a=c(-1,1)) {
e1=1.0-0.5*(e/100.0)
Z=qnorm(e1)^2
aux=FindLambda.MA(x,n,Z=Z,a=a); l1ze0=aux[1]; l2ze0=aux[2]
aux=FindLambda.MA(x,n,Z=Z,a=a,correct=T); l1ze0c=aux[1]; l2ze0c=aux[2]
pp=x/n
I=c(0,0);S=c(0,0)
if(pp[1]==0) {I[1]=1} else {I[1]=0}
if(pp[2]==1) {I[2]=1} else {I[2]=0}
if(pp[1]==1) {S[1]=1} else {S[1]=0}
if(pp[2]==0) {S[2]=1} else {S[2]=0}
Z=qnorm(e1)
h=Z^2*(1.0+2.0*I)/4.0
xx=x+h; nn=n+2*h; pp=xx/nn
fact=sum(a^2*(pp*(1-pp))/nn) ; l1zw4=a%*%pp -Z*sqrt(fact) # lower limit
h=Z^2*(1.0+2.0*S)/4.0
xx=x+h; nn=n+2*h; pp=xx/nn
fact=sum(a^2*(pp*(1-pp))/nn); l2zw4=a%*%pp +Z*sqrt(fact)
xx=x+0.5; nn=n+1; b=sqrt(xx/nn)
fact1=asin(b[2])-asin(b[1])
valor=abs(a[2])
aux=uniroot.all(As, c(-1, 1),x=xx,n2=nn,a=a,fact1=fact1,e1=e1,maxiter = 2000, n = 5000)
lim_inf1=min(aux)*valor; lim_sup1=max(aux)*valor
aux=uniroot.all(AM, c(-1, 1),x=xx,n2=nn,a=a,e1=e1, maxiter = 2000, n = 5000)
lim_inf2=min(aux)*valor; lim_sup2=max(aux)*valor
# xx=x+Z^2/4; nn=n+Z^2/2
# aux=FindLambda.MA(x=xx,n=nn,Z=Z^2,correct=F,a);
# lim_inf3=min(aux)*valor; lim_sup3=max(aux)*valor
# aux=FindLambda.MA(x=xx,n=nn,Z=Z^2,correct=T,a);
# lim_inf4=min(aux)*valor; lim_sup4=max(aux)*valor
# l1ze3=lim_inf3;l2ze3=lim_sup3;
# l1ze3c=lim_inf4;l2ze3c=lim_sup4;
# l1ze3=max(-1.0,l1ze3);l2ze3=min(1.0,l2ze3)
# l1ze3c=max(-1.0,l1ze3c);l2ze3c=min(1.0,l2ze3c)
l1ae1=lim_inf1;l2ae1=lim_sup1;
l1aem1=lim_inf2;l2aem1=lim_sup2;
l1ze0=max(-1.0,l1ze0);l2ze0=min(1.0,l2ze0)
l1ze0c=max(-1.0,l1ze0c);l2ze0c=min(1.0,l2ze0c)
l1ae1=max(-1.0,l1ae1);l2ae1=min(1.0,l2ae1)
l1aem1=max(-1.0,l1aem1);l2aem1=min(1.0,l2aem1)
l1zw4=max(-1.0,l1zw4);l2zw4=min(1.0,l2zw4)
result=cbind(
c(l1ze0,l1ze0c,l1zw4,l1ae1,l1aem1),
c(l2ze0,l2ze0c,l2zw4,l2ae1,l2aem1))
#l1ze3,l1ze3c,
#l2ze3,l2ze3c,
rownames(result)=c("Score (without cc)","Score (with cc)","Adjusted Wald","Adjusted Arc Sine","Adjusted-M Arc Sine")
colnames(result)=c("lower limit", "upper limit")
# "Adjusted Score (without cc)", "Adjusted Score (with cc)",
result}
test2.MA=function(x,n,a,lambda=0,e=5,
alternative = c("two.sided", "less", "greater")) {
alternative <- match.arg(alternative)
ze0=FindZ.MA(x,n,correct=FALSE,a=a,lambda=lambda,plot=TRUE)
ze0c=FindZ.MA(x,n,correct=TRUE,a=a,lambda=lambda,plot=TRUE)
e1=1.0-(e/100.0)/2.0;Z=qnorm(e1)
xx=x+Z^2/4;nn=n+Z^2/2
p=xx/nn; d=a%*%p
#pp=sum(xx)/sum(nn);qq=1.0-pp;
pp=Pvero_dif(delta=lambda/a[2],x=xx,n=nn); qq=1-pp
fact1=a[1]^2*pp*qq/nn[1]+a[2]^2*(pp+lambda)*(1-pp-lambda)/nn[2]; fact2=nn[1]*nn[2]+nn[1]+nn[2];
if (nn[1]==nn[2]) {c=2.0/fact2} else {c=1.0/fact2}
ze3=sqrt((d-lambda)^2/fact1) # any value of lambda (not only 0)
if (d-lambda>0) {dc=d-lambda} else {dc=-d+lambda}
if (dc>c) {e3c=(dc-c)^2/fact1} else {e3c=0.0}
ze3c=sqrt(e3c)
ze0=sqrt(ze0);ze0c=sqrt(ze0c)
pp=x/n
I=c(0,0);S=c(0,0)
if(pp[1]==0) {I[1]=1} else {I[1]=0}
if(pp[2]==1) {I[2]=1} else {I[2]=0}
if(pp[1]==1) {S[1]=1} else {S[1]=0}
if(pp[2]==0) {S[2]=1} else {S[2]=0}
dif=a%*%pp
if(dif>lambda)
{h=Z^2*(1.0+2.0*I)/4.0} else
{h=Z^2*(1.0+2.0*S)/4.0}
xx=x+h; nn=n+2*h; pp=xx/nn
fact=sum(a^2*(pp*(1-pp))/nn)
zw4=(a%*%pp -lambda)/sqrt(fact)
xx=x+0.5; nn=n+1; b=sqrt(xx/nn)
fact1=asin(b[2])-asin(b[1])
# aux=uniroot.all(AM, c(-1, 1),x=xx,n2=nn,a=a,e1=e1,maxiter = 2000, n = 5000)
# no necesary
p_vero1=Pvero_dif(delta=lambda/a[2],x=xx,n=nn)
p_vero2=p_vero1+lambda/a[2]
if(p_vero2<0.0) p_vero2=0.0; if(p_vero2>1.0) p_vero2=1.0
b=c(b,sqrt(p_vero1),sqrt(p_vero2))
fact2=asin(b[2])-asin(b[1])-asin(b[4])+asin(b[3]);
ae1=4.0*nn[1]*nn[2]*fact2^2/(nn[1]+nn[2])
ae1m=4.0*(nn[1]* ( asin(b[1]) - asin(b[3]) )^2 + nn[2]* (asin(b[2]) - asin(b[4]) ) ^2)
cc3=sqrt(ae1);cc5=sqrt(ae1m);cc4=zw4;
zs=c( ze0, ze0c, cc4, cc3, cc5);zs=abs(zs)
# ze3,ze3c,
pvalue=numeric(5)
if(a%*%(x/n)<lambda){zs=(-zs)} # eran todos positivos
pvalue=switch(alternative,two.sided = 2.0*pnorm(-abs(zs)),less = pnorm(zs),greater = 1-pnorm(zs))
result=cbind(zs,pvalue)
rownames(result)=c("Score (without cc)","Score (with cc)","Adjusted Wald","Adjusted Arc Sine", "Adjusted-M Arc Sine")
colnames(result)=c("z value", "p-value")
# "Adjusted Score (without cc)", "Adjusted Score (with cc)",
result}
prop.test2.MA=function(x,n,lambda,conf.level=0.95,
alternative = c("two.sided", "less", "greater")) {
a=c(-1,1)
e=100*(1-conf.level)
A=test2.MA(x,n,a,lambda,e,alternative)
alternative <- match.arg(alternative)
B=switch(alternative,two.sided = ic2.MA(x,n,e,a),less = ic2.MA(x,n,e=2*e,a),greater = ic2.MA(x,n,e=2*e,a))
B[,1]=switch(alternative,two.sided = B[,1],less = -1,greater=B[,1])
B[,2]=switch(alternative,two.sided = B[,2],less = B[,2],greater=1)
result=cbind(B,A)
list(estimate=x/n,a=a,difference=as.numeric(a%*%(x/n)),inference=result,alternative=alternative,lambda=lambda)
}
#### Combination of two proportions: testComb2.MA + icComb2.MA
###############################
testComb2.MA=function(x,n,a,lambda=0,e=5,
alternative = c("two.sided", "less", "greater")) {
alternative = match.arg(alternative)
prop=x/n; L=a%*%prop
ze0=FindZ.MA(x,n,correct=FALSE,a=a,lambda=lambda)
ze0c=FindZ.MA(x,n,correct=TRUE,a=a,lambda=lambda)
xx=x+1;nn=n+2; p=xx/nn; L=a%*%p; d=L-lambda
fact=sum(a^2*p*(1-p)/nn); zw2=d^2/fact;
e1=1.0-0.5*(e/100);Z=qnorm(e1)
if (a[1]>0) {s1=1.0} else {s1=-1.0}
if (a[2]>0) {s2=1.0} else {s2=-1.0}
pp11=0.5*(1.0+s1);pp12=0.5*(1.0+s2)
pp21=0.5*(1.0-s1);pp22=0.5*(1.0-s2)
if (pp11==prop[1]) {I1=1} else {I1=0}
if (pp12==prop[2]) {I2=1} else {I2=0}
if (pp21==prop[1]) {S1=1} else {S1=0}
if (pp22==prop[2]) {S2=1} else {S2=0}
if (a%*%prop>lambda)
{h=c(1+I1*2,1+I2*2)} else {h=c(1+S1*2,1+S2*2)}
xxx=x+Z^2*h/4; nnn=n+Z^2*h/2; pd=xxx/nnn; Ld=a%*%pd; dd=Ld-lambda
factd=sum(a^2*pd*(1-pd)/nnn)
zw4=dd^2/factd
# case LR1*/
xm=x+0.5; nm=n+1;pm=xm/nm
#Calculation of p_verosimilitud with increase data*/
ad1=sum(xm)
alpha=lambda/a[2];bet=-a[1]/a[2]
c0=alpha*(1.0-alpha)*xm[1]
c1=bet*ad1-nm[1]*alpha*(1.0-alpha)-alpha*bet*(nm[2]+2.0*xm[1])
c2=bet*(alpha*(nm[2]+2*nm[1])-(nm[1]+xm[2])-bet*(nm[2]+xm[1]))
c3=bet*bet*(nm[1]+nm[2])
A=-4.5*c3*(c1*c2-3.0*c0*c3)+c2*c2*c2
B=c2*c2-3.0*c1*c3
if (A>=0) {u=sqrt(B)} else {u=-sqrt(B)}
w1= A/(u^3); if (w1>1.0) w1=1.0; if (w1<(-1.0)) w1=-1.0
fhi=(pi+ acos(w1))/3
pvm1=(-1*c2+2*cos(fhi)*u)/(3*c3) ;if (pvm1>1) pvm1=1; if (pvm1<0) pvm1=0
pvm2=alpha+bet*pvm1;if (pvm2>1) pvm2=1;if (pvm2<0) pvm2=0;
gk1=log(pvm1/pm[1]);gk2=log((1-pvm1)/(1-pm[1]));gk3=log(pvm2/pm[2]);gk4=log((1-pvm2)/(1-pm[2]));
lr1=-2*(xm[1]*gk1+(nm[1]-xm[1])*gk2+xm[2]*gk3+(nm[2]-xm[2])*gk4)
cc=c(sqrt(ze0),sqrt(ze0c),sqrt(zw2),sqrt(zw4),sqrt(lr1))
if(a%*%(x/n)<lambda){cc=-cc} # eran todos positivos
pvalor=switch(alternative,two.sided = 2.0*pnorm(-abs(cc)),less = pnorm(cc),greater = 1-pnorm(cc))
result=cbind(cc,pvalor)
rownames(result)=c("Score (without cc)","Score (with cc)",
"Adjusted Wald(a)", "Adjusted Wald(b)" ,"Adjusted likelihood ratio" )
colnames(result)=c("z value", "p-value")
result
}
icComb2.MA=function(x,n,e=5,a) {
prop=x/n; e1=1-(e/100)/2; Z=qnorm(e1)^2
limites=c(sum(a[a<0]),sum(a[a>0]))
lze0=FindLambda.MA(x,n,Z=Z,correct=FALSE,a=a)
lze0c=FindLambda.MA(x,n,Z=Z,correct=TRUE,a=a)
xx=x+1;nn=n+2; p=xx/nn; L=a%*%p; fact=sum(a^2*p*(1-p)/nn)
Z=qnorm(e1)
lzw2=c(L-Z*sqrt(fact),L+Z*sqrt(fact))
if (a[1]>0) {s1=1.0} else {s1=-1.0}
if (a[2]>0) {s2=1.0} else {s2=-1.0}
pp11=0.5*(1.0+s1);pp12=0.5*(1.0+s2)
pp21=0.5*(1.0-s1);pp22=0.5*(1.0-s2)
if (pp11==prop[1]) {I1=1} else {I1=0}
if (pp12==prop[2]) {I2=1} else {I2=0}
if (pp21==prop[1]) {S1=1} else {S1=0}
if (pp22==prop[2]) {S2=1} else {S2=0}
#limite inferior
h=c(1+I1*2,1+I2*2)
xxx=x+Z^2*h/4; nnn=n+Z^2*h/2; pd=xxx/nnn; Ldi=a%*%pd
factdi=sum(a^2*pd*(1-pd)/nnn)
#limite superior
h=c(1+S1*2,1+S2*2)
xxx=x+Z^2*h/4; nnn=n+Z^2*h/2; pd=xxx/nnn; Lds=a%*%pd
factds=sum(a^2*pd*(1-pd)/nnn)
lzw4=c(Ldi-Z*sqrt(factdi),Lds+Z*sqrt(factds))
# /*caso de LR1*/
xm=x+0.5;nm=n+1;pm=xm/nm
llr=uniroot.all(RazonVero, limites,x=xm,n2=nm,a=a,Z=Z,maxiter = 2000, n = 5000)
result=rbind(lze0,lze0c,lzw2,lzw4,llr)
rownames(result)=c("Score (without cc)","Score (with cc)","Adjusted Wald (a)",
"Adjusted Wald (b)", "Adjusted likelihood ratio" )
colnames(result)=c("lower limit", "upper limit")
result[result[,1]<limites[1]]=limites[1]; result[result[,1]>limites[2]]=limites[2]
result}
prop.Comb2.MA=function(x,n,a,lambda,conf.level=0.95,
alternative = c("two.sided", "less", "greater")) {
e=100*(1-conf.level)
A=testComb2.MA(x,n,a,lambda,e,alternative)
alternative =match.arg(alternative)
B=switch(alternative,two.sided = icComb2.MA(x,n,e,a),less = icComb2.MA(x,n,e=2*e,a),
greater = icComb2.MA(x,n,e=2*e,a))
B[,1]=switch(alternative,two.sided = B[,1],less = sum(a[a<0]),greater=B[,1])
B[,2]=switch(alternative,two.sided = B[,2],less = B[,2],greater=sum(a[a>0]))
result=cbind(B,A)
list(estimate=x/n,a=a,L=as.numeric(a%*%(x/n)),inference=result,alternative=alternative,lambda=lambda)
}
#### Combination of many proportions: testCombi.MA + icCombi.MA
###############################
testCombi.MA=function(x,n,a,lambda=0,e=5,
alternative = c("two.sided", "less", "greater")) {
alternative = match.arg(alternative)
ze0=FindZ.MA(x,n,correct=FALSE,a=a,lambda=lambda,plot=TRUE)
ze0c=FindZ.MA(x,n,correct=TRUE,a=a,lambda=lambda,plot=TRUE)
e1=1.0-0.5*(e/100); Z=qnorm(e1)
fact1=0.5*Z^2; pr=x/n ;L=a%*%pr
K=length(pr); prop_sa=numeric(K)
s=rep(1,K);s[a<0]=-1
prop_sa=0.5*(1+s); prop_sb=0.5*(1-s)
h=numeric(K)
if (L>lambda){ii=prop_sa==pr; h[ii]=fact1*(1+K)/K; h[!ii]=fact1/K}
else {ii=prop_sb==pr; h[ii]=fact1*(1+K)/K; h[!ii]=fact1/K}
xx=x+h;nn=n+2*h;pp=xx/nn; L=a%*%pp
denom=sum(a^2*pp*(1-pp)/nn); denom2=sum(a^2/n)
zw3=sqrt((L-lambda)^2/denom)
B=sum(a); denom3=denom2-(B-2*lambda)^2/sum(n)
pa0=4.0*(a%*%pr-lambda)^2/denom3; zpa0=sqrt(pa0)
z=numeric(4)
z[1]=sqrt(ze0); z[2]=sqrt(ze0c);z[3]=zw3;z[4]=zpa0
if(a%*%(x/n)<lambda){z=-z} # eran todos positivos
pvalue=switch(alternative,two.sided = 2.0*pnorm(-abs(z)),less = pnorm(z),greater = 1-pnorm(z))
result=cbind(z,pvalue)
rownames(result)=c("Score (without cc)","Score (with cc)",
"Adjusted Wald", "Peskun "); colnames(result)=c("z value", "p-value")
result}
icCombi.MA=function(x,n,e=5,a) {
e1=1.0-0.5*(e/100.0);Z=qnorm(e1)^2
aux=FindLambda.MA(x,n,Z=Z,a=a); l1ze0=aux[1]; l2ze0=aux[2]
aux=FindLambda.MA(x,n,Z=Z,correct=T,a=a); l1ze0c=aux[1]; l2ze0c=aux[2]
fact1=0.5*Z; pr=x/n ;L=a%*%pr
K=length(pr); prop_sa=numeric(K)
s=rep(1,K);s[a<0]=-1
prop_sa=0.5*(1+s); prop_sb=0.5*(1-s)
h1=numeric(K);h2=numeric(K)
ii=prop_sa==pr; h1[ii]=fact1*(1+K)/K; h1[!ii]=fact1/K
ii=prop_sb==pr; h2[ii]=fact1*(1+K)/K; h2[!ii]=fact1/K
xx1=x+h1;nn1=n+2*h1;pp1=xx1/nn1; L1=a%*%pp1
xx2=x+h2;nn2=n+2*h2;pp2=xx2/nn2; L2=a%*%pp2
denom1=sum(a^2*pp1*(1-pp1)/nn1)
denom2=sum(a^2*pp2*(1-pp2)/nn2)
lim=matrix(ncol=2,nrow=4)
lim[1,]=c(l1ze0,l2ze0); lim[2,]=c(l1ze0c,l2ze0c)
lim[3,]=c(L1-qnorm(e1)*sqrt(denom1),L2+qnorm(e1)*sqrt(denom2))
B=sum(a);sumn=sum(n); c1=sumn/(sumn+Z);
L2=a%*%pr; c2=L2+B*Z/(2*sumn);c3=(B-2*L2)^2/sumn;
c4=sum(a^2/n)*(sumn+Z)/sumn;c5=c4-c3;
lim[4,]=c(c1*(c2-qnorm(e1)*sqrt(c5)*0.5),c1*(c2+qnorm(e1)*sqrt(c5)*0.5))
B_neg=sum(a[a<0]);B_pos=sum(a[a>0])
ii=lim[,1]<B_neg; lim[ii,1]=B_neg
ii=lim[,2]>B_pos; lim[ii,2]=B_pos
rownames(lim)=c("Score (without cc)","Score (with cc)","Adjusted Wald",
"Peskun ")
colnames(lim)=c("lower limit", "upper limit")
lim}
prop.Combi.MA=function(x,n,a,lambda,conf.level=0.95,
alternative = c("two.sided", "less", "greater")) {
e=100*(1-conf.level)
A=testCombi.MA(x,n,a,lambda,e,alternative)
alternative =match.arg(alternative)
B=switch(alternative,two.sided = icCombi.MA(x,n,e,a),less = icCombi.MA(x,n,e=2*e,a),
greater = icCombi.MA(x,n,e=2*e,a))
B[,1]=switch(alternative,two.sided = B[,1],less = sum(a[a<0]),greater=B[,1])
B[,2]=switch(alternative,two.sided = B[,2],less = B[,2],greater=sum(a[a>0]))
result=cbind(B,A)
list(estimate=x/n,a=a,L=as.numeric(a%*%(x/n)),inference=result,alternative=alternative,
lambda=lambda)
}
###########################
Pvero_dif=function(delta,x,n){
a1=sum(x); a=sum(n)
b=-(a+a1-(a+n[1])*delta)
c=-delta*(-n[1]*delta+2.0*x[1]+sum(n))+a1
d=x[1]*delta*(1-delta)
v1=b/(3.0*a);v2=v1^3;v3=(b*c)/(6.0*a^2);v4= d/(2.0*a)
v=v2-v3+v4
u= sqrt( b^2/(9.0*a^2)-c/(3*a))
ii=v<0; u[ii]=-u[ii]
w1= v/(u^3)
w1[w1>1]=1.0; w1[w1<(-1.0)]=-1.0
w2=acos(w1); w=(acos(-1.0)+w2)/3.0
pvero=2.0*u*cos(w)-b/(3*a)
pvero[pvero>1.0]=1.0; pvero[pvero<0.0]=0
pvero}
###########################
# END OF AUXILIAR FUNCTIONS
############################
OK=complete.cases(x, n); x=x[OK]; n=n[OK]
k=length(x)
recomen=""
##############################################################
# A PROPORTION
if (k==1) {
if (missing(p)) p=0.5
else if (p<=0 | p>=1) stop("value of 'p' must be in (0,1)")
resultado=prop.test.MA(x,n,p0=p,conf.level,alternative)
# Recommendation ____
if (conf.level==0.90 & sum(n)>=100) {recomen="Adjusted Arc Sine"}
else if (conf.level!=0.90 & sum(n)<=80) {recomen="Modified Score"}
else {recomen= "Adjusted Wald"} }
##############################################################
# DIFFERENCE OF PROPORTIONS
else if(k==2 & missing(a)) {
if (missing(p)) p=0
else if (p<=-1 | p>=1) stop("value of 'p' must be in (-1,1)")
resultado=prop.test2.MA(x=x,n=n,lambda=p,conf.level,alternative)
a=resultado$a
# Recommendation ____
#if (p==0) {
# if (conf.level==0.90) recomen="Adjusted Score (without cc)"
# if (conf.level!=0.90 & n[1]==n[2]) recomen="Adjusted Score (with cc)"
#if (conf.level!=0.90 & n[1]!=n[2]) recomen="Score (with cc)"
#else recomen="No recommendation" }
#else {
if (conf.level==0.99) recomen="Adjusted-M Arc Sine"
if (conf.level!=0.99 & sum(n)>=200 & n[1]==n[2]) recomen="Adjusted Wald"
if (conf.level!=0.99 & (sum(n)<200 | n[1]!=n[2])) recomen="Adjusted Arc Sine"
#}
}
##############################################################
# COMBINATION OF TWO PROPORTIONS
else if(k==2 & !missing(a)) {
if (missing(p)) p=0.5*(sum(a[a<0])+ sum(a[a>0]) )
else {
ai=sum(a[a<0]); as=sum(a[a>0])
if (p<=ai | p>=as) {
cat(paste("a_neg=",ai,": sum of negative coefficients",sep=""),"\n")
cat(paste("a_pos=",as,": sum of positive coefficients",sep=""),"\n")
stop("value of 'p' must be in (a_neg,a_pos)")
}}
resultado=prop.Comb2.MA(x=x,n=n,a=a,lambda=p,conf.level,alternative)
# Recommendation ____
if (conf.level!=0.90 & n[1]<60 & n[2]<60) recomen="Adjusted Wald(a)"
if (conf.level!=0.99) recomen="Adjusted Wald(b)"
if (conf.level==0.99) recomen="Adjusted Likelihood Ratio"
}
##############################################################
# COMBINATION OF MORE THAN TWO PROPORTIONS
else if(k>2) {
if (missing(a)) stop ("value of 'a' is lost")
if (missing(p)) p=0.5*(sum(a[a<0])+ sum(a[a>0]) )
else {
ai=sum(a[a<0]); as=sum(a[a>0])
if (p<=ai | p>=as) {
cat(paste("a_neg=",ai,": sum of negative coefficients",sep=""),"\n")
cat(paste("a_pos=",as,": sum of positive coefficients",sep=""),"\n")
stop("value of 'p' must be in (a_neg,a_pos)")}}
resultado=prop.Combi.MA(x=x,n=n,a=a,lambda=p,conf.level,alternative)
if (k==3 & max(n)>10) recomen="Score (with cc)"
else if (k>3 & max(n)>10) recomen="Score (without cc)"
else if ( max(n)<=10) recomen="Adjusted Wald"
}
if (conf.level!=0.90 & conf.level!=0.95 &conf.level!=0.99) {recomen="No recommendation"}
#if (recomen=="") recomen="No recommendation"
resultado$k=k; resultado$x=x; resultado$n=n;
resultado$a=a; resultado$p=p; resultado$conf.level=conf.level;
resultado$recommendation=recomen
#####################################################
# noa=missing(a)
# a0=a
if (coverage)
{
if (noa) {resultado=coveragecomb(resultado,nrep)}
else {resultado=coveragecomb(resultado,nrep,a=a0)}
}
object=resultado
class(object) ="prop.comb"
object
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.