R help - Dec 2012 - R function for computing Simultaneous confidence intervals for multinomial proportions

Dear all,

Does someone know an R function implementing the method of Sison and 
Glaz (1995) (see full ref below) for computing Simultaneous confidence 
intervals for multinomial proportions?

As alternative method, I think to boostrap the mean of each proportion 
and get in that way confidence interval of the mean.

I observed 21 times a response that could be one out of 8 categories 
(multinomial response). I computed the proportions for every categories.
I did it independently 12 times. Hence I have 12 replicats for each 
proportions.
Is boostraping the mean proportion over the 12 replicats for getting the 
confidence interval is a good idea?

I tried (see codes below) and gets some confidence interval quiet large. 
Actually, according to the boostraped CI, there are no difference in 
proportions, while a chi2 test confirms that there are.

############################################ generation of a data set: 
multinomial réponse variable repeatedly measured
sub<-c("s1","s2","s3","s4")#subject
per<-c("p1","P2","p3")#period, RM within
subject
tim<-paste("t",1:21,sep="")#tim, RM within period
mydata<-expand.grid(tim=tim,per=per,sub=sub)
#for simplicity, let's consider that period is not a repeated measure 
within subject. Hence, we have 12 trials ("independent") of 21
observations

#random generation of the response variable
require(plyr)
posl<-paste("z",1:8,sep="")
mydata$pos<-factor(NA,levels=posl)
n=nrow(mydata)
mydata$pos<-replicate(n,sample(posl,1,prob=c(0.05,0.05,0.05,0.4,0.3,0.05,0.05,0.05)))#strong
preference for z3 and z4

##pre-treatment of the above row data, number of time (out of 21) each 
position was observed
mydata2=ddply(mydata,.(sub,per,pos),summarize,nobs=length(pos),.drop=F)
mydata2[,1:3]=catcolwise(function(x)as.factor(x))(mydata2)
summary(mydata2)

# boxplot of frequencies of occpupancy
require(ggplot2)
mydata2$fobs=mydata2$nobs/21
ggplot(mydata2)+geom_boxplot(aes(pos,fobs))

###chi2 test
nobsT=ddply(mydata,.(pos),summarize,T=length(pos))
nobsT=nobsT[,2]
chisq.test(nobsT)

#######################bootstraping of the mean proportions over the 12 
replicats
require(boot)
maf=function(data,index){o=ddply(data[index,],.(pos),summarize,mmea=mean(fobs),.drop=F);as.vector(o[,2])}
bootobj=boot(mydata2,maf,R=9999,stype="i")

confint=ddply(mydata2,.(pos),summarize,mean.prop=mean(fobs))
confint$boot.LOW=-999
confint$boot.UP=-999
for (posi in 1:8){
try({confint$boot.LOW[posi]<-boot.ci(bootobj,0.95,index=posi,type="bca")$bca[4]},silent=T)
try({confint$boot.UP[posi]<-boot.ci(bootobj,0.95,index=posi,type="bca")$bca[5]},silent=T)
}

####################### plotting boostraped CI
ggplot(mydata2)+geom_boxplot(aes(pos,fobs))+
geom_errorbar(data=confint,aes(pos,ymin=mean.prop-boot.LOW,ymax=mean.prop+boot.UP),lty=2,width=0.5)+
   geom_point(data=confint,aes(pos,mean.prop),shape=24,size=2)


References of the CI estimations I would like tu run:

Glaz, J. and C.P. Sison. 1999. Simultaneous confidence intervals for 
multinomial proportions. Journal of Statistical Planning and Inference 
82:251-262.

Sison, C.P and J. Glaz. 1995. Simultaneous confidence intervals and 
sample size determination for multinomial proportions. Journal of the 
American Statistical Association, 90:366-369. Paper available at 
http://tx.liberal.ntu.edu.tw/~purplewoo/Literature/!Methodology/!Distribution_SampleSize/SimultConfidIntervJASA.pdf

I thank you in advance,
Respectfully,
-- 
*Pierre THIRIET*

Doctorant en écologie marine - Marine Ecology PhD student
----------------------------------------------------------------------------------------------------------------------------------------
EA 4228 - ECOMERS - /Ecosystèmes Côtiers Marins et Réponses aux Stress/
Université de Nice - Sophia Antipolis, Faculté des Sciences, Parc 
Valrose, 06108 Nice Cedex 2, France
More about my work and my lab. 
<http://www.unice.fr/ecomers/index.php?option=com_comprofiler&task=userProfile&user=101&Itemid=111&lang=fr>
----------------------------------------------------------------------------------------------------------------------------------------
Cell: +336 79 44 91 90
Office: +334 92 07 68 33
Skype: pierre.d.thiriet
Mail: pierre.thiriet@unice.fr
Mail(bis): pierre.d.thiriet@gmail.com

	[[alternative HTML version deleted]]