simulaAnova=function(n=100,ngrupos=2,efecto=10,media=170,sd=10){ x=rnorm(n*ngrupos,media,sd); grupos=as.factor(rep(c("A","B","C","D","E")[1:ngrupos],n)); x[grupos=="A"]=x[grupos=="A"]+efecto; beanplot(x ~ as.factor(grupos),cut=1,col="brown",ll=0.04,ylim=c(media-3*sd,media+3*sd)) } simulaChi2=function(ncol=c(100,30),p=c(.5,.5)){ #ncol=c(100,30) #p=c(.5,.5) v1=c() for(i in 1:length(ncol)){v1=c(v1,rep(i,ncol[i]))} v2=v1 for(i in 1:length(v1)){ v2[i]=1+rbinom(1,1,p[v1[i]]) } v1=as.factor(v1) v2=as.factor(v2) levels(v1)=paste("Group",c("A","B","C","D","E","F","G","H")[1:length(ncol)],sep=" ") levels(v2)=c("No","Yes") tabla=table(v1,v2) mosaicplot(tabla,col=c("lightblue","royalblue"),sub=sprintf("p=%1.3f",chisq.test(tabla)$p.value),ylab="Something interesting happens?",xlab="",off=.4) resultados=round(tabla[,2]/(apply(tabla,1,sum)),3) names(resultados)=levels(v1) xcoord=rep(0,length(ncol)) for(i in 2:length(ncol)){xcoord[i]=xcoord[i-1]+ncol[i-1]} xcoord=xcoord/sum(ncol) #print(xcoord) print(resultados) text(xcoord,rep(0.005,length(ncol)),paste(100*resultados,"%",sep=""),cex=1.2,col="yellow",pos=4) } simulaChi2(c(150,50)) #################### simulaChi2(c(150,25,25),c(.5,.5,.5)) ################ library(foreign) pm=read.spss("penamuerte-chi2-confusion.sav",T,T) pm$Victima_Acusado=sprintf("%s-%s",pm$rzVictima,pm$rzAcusado) pm$Acusado_Victima=sprintf("%s-%s",pm$rzAcusado,pm$rzVictima) tabla=xtabs( ~ rzAcusado +penaM,data=pm) mosaicplot(tabla,col=c("lightblue","blue"),shade=F,off=1,las=2) tabla=xtabs( ~ rzVictima +penaM,data=pm) mosaicplot(tabla,col=c("lightblue","blue"),shade=F,off=1,las=2) tabla=xtabs( ~ rzAcusado +rzVictima,data=pm) mosaicplot(tabla,col=c("brown","DarkSalmon "),shade=F,off=1,las=2) #tabla=xtabs( ~ rzAcusado +penaM+rzVictima,data=pm) tabla=xtabs( ~ Acusado_Victima+penaM,data=pm) mosaicplot(tabla,col=c("lightblue","blue"),shade=F,off=1,las=2) tabla=xtabs( ~ Victima_Acusado+penaM,data=pm) mosaicplot(tabla,col=c("lightblue","blue"),shade=F,off=1,las=2) ########## library(psych) extraeXY=function(nind=1000,coef=.9,ny=2){ x=seq(0,1,length.out=nind) y1=x+coef*rnorm(nind,0,1-coef) res=data.frame(x,y1) for(i in 2:ny){ yi=x+rnorm(nind,0,(1-coef)) res=cbind(res,yi) } names(res)=c("x",sprintf("y%d",1:ny)) res } alpha(extraeXY(ny=10,coef=(.2),nind=200)[,-1])