R/twobici.R

twobici <-
function(r1=sum(x),n1=length(x),r2=sum(y),n2=length(y),x=NA,y=NA,alpha=.05){
#
# Compute confidence interval for p1-p2,
# the difference between probabilities of
# success for a two binomials using Beal's method.
#
# r is number of successes
# n is sample size
# if x contains data, r1 is taken to be the
# number of 1s in x and n1 is length(x)
#
if(length(r1)>1)stop("r1 must be a single number, not a vector")
if(length(n1)>1)stop("n1 must be a single number, not a vector")
if(length(r2)>1)stop("r2 must be a single number, not a vector")
if(!is.na(sum(r1)) || !is.na(sum(n1)) || !is.na(sum(r2)) || !is.na(sum(n2))){
if(r1<0 || n1<0)stop("Both r1 and n1 must be greater than 0")
if(r1 > n1)stop("r1 can't be greater than n1")
if(r2<0 || n2<0)stop("Both r2 and n2 must be greater than 0")
if(r2 > n2)stop("r2 can't be greater than n2")
}
if(!is.na(sum(x))){
r1<-sum(x)
n1<-length(x)
}
if(!is.na(sum(y))){
r2<-sum(y)
n2<-length(y)
}
a<-(r1/n1)+(r2/n2)
b<-(r1/n1)-(r2/n2)
u<-.25*((1/n1)+(1/n2))
v<-.25*((1/n1)-(1/n2))
V<-u*((2-a)*a-b^2)+2*v*(1-a)*b
crit<-qchisq(1-alpha/2,1)
A<-sqrt(crit*(V+crit*u^2*(2-a)*a+crit*v^2*(1-a)^2))
B<-(b+crit*v*(1-a))/(1+crit*u)
ci<-NA
ci[1]<-B-A/(1+crit*u)
ci[2]<-B+A/(1+crit*u)
p1<-r1/n1
p2<-r2/n2
list(ci=ci,p1=p1,p2=p2)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.