# dat <- read.delim("http://phi.med.gunma-u.ac.jp/msb/data/p11.txt") # attach(dat) # X <- table(GRP,PARITY) # HF <- X[1,] # MF <- X[2,] # PF <- X[3,] # 以上のコードは,ファイルからデータを読まずに作るには以下で代替可 HF <- c(4,7,7,6,4,1,0,1,0,0) names(HF) <- 0:9 H <- rep(0:9,HF) MF <- c(1,3,6,6,9,3,0,1,1,0) names(MF) <- 0:9 M <- rep(0:9,MF) PF <- c(1,2,4,8,7,2,3,1,1,1) names(PF) <- 0:9 P <- rep(0:9,PF) PARITY <- c(H,M,P) GRP <- factor(c(rep('H',30),rep('M',30),rep('P',30))) dat <- data.frame(PARITY=PARITY,GRP=GRP) # ここまで win.metafile("it11-ans-2006-1.emf",width=6,height=6,pointsize=14) par(family="sans",mai=c(0.8,0.5,0.4,0.5),las=1) boxplot(PARITY~GRP,main="地域別生涯出生児数") dev.off() fligner.test(PARITY~GRP) kruskal.test(PARITY~GRP) pairwise.wilcox.test(PARITY,GRP,exact=F) # 以上が普通の回答 library(vcd) win.metafile("it11-ans-2006-2.emf",width=6.7,height=10,pointsize=14) par(family="sans",mai=c(0.4,0.4,0.4,0.4),las=1,mfrow=c(3,1)) XHP <- goodfit(H,"poisson") SXHP <- summary(XHP) TXHP <- paste("ポアソン分布 (p=",sprintf("%4.2f",SXHP[3]),")") XHN <- goodfit(H,"nbinom") SXHN <- summary(XHN) TXHN <- paste("負の2項分布 (p=",sprintf("%4.2f",SXHN[3]),")") ix <- barplot(HF,main="H市の出生児数分布と分布の当てはめ",ylim=c(0,10)) lines(ix,predict(XHP,newcount=0:9),lty=1,col="red") lines(ix,predict(XHN,newcount=0:9),lty=2,col="blue") legend(8,max(HF),lty=c(1,2),legend=c(TXHP,TXHN),col=c("red","blue")) XMP <- goodfit(M,"poisson") SXMP <- summary(XMP) TXMP <- paste("ポアソン分布 (p=",sprintf("%4.2f",SXMP[3]),")") XMN <- goodfit(M,"nbinom") SXMN <- summary(XMN) TXMN <- paste("負の2項分布 (p=",sprintf("%4.2f",SXMN[3]),")") ix <- barplot(MF,main="M村の出生児数分布と分布の当てはめ",ylim=c(0,10)) lines(ix,predict(XMP,newcount=0:9),lty=1,col="red") lines(ix,predict(XMN,newcount=0:9),lty=2,col="blue") legend(8,max(MF),lty=c(1,2),legend=c(TXMP,TXMN),col=c("red","blue")) XPP <- goodfit(P,"poisson") SXPP <- summary(XPP) TXPP <- paste("ポアソン分布 (p=",sprintf("%4.2f",SXPP[3]),")") XPN <- goodfit(P,"nbinom") SXPN <- summary(XPN) TXPN <- paste("負の2項分布 (p=",sprintf("%4.2f",SXPN[3]),")") ix <- barplot(PF,main="P村の出生児数分布と分布の当てはめ",ylim=c(0,10)) lines(ix,predict(XPP,newcount=0:9),lty=1,col="red") lines(ix,predict(XPN,newcount=0:9),lty=2,col="blue") legend(8,max(PF),lty=c(1,2),legend=c(TXPP,TXPN),col=c("red","blue")) dev.off()