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.thir...@unice.fr
Mail(bis): pierre.d.thir...@gmail.com
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
