dat=read.csv("mendelPeaCp.csv", as.is = T, header = T)
layout(matrix(c(1,2,3), 1, 3, byrow = T))
dat$Realobs1=dat$Total*(dat$MendelRatio1/(dat$MendelRatio1+dat$MendelRatio2))
dat$Realobs2=dat$Total*(dat$MendelRatio2/(dat$MendelRatio1+dat$MendelRatio2))
dat$MendelChi=(dat$observed1-dat$Realobs1)^2/dat$Realobs1+(dat$observed2-dat$Realobs2)^2/dat$Realobs2
od=order(dat$MendelChi)
xobs=qchisq(seq(1/nrow(dat), 1-1/nrow(dat), length.out = nrow(dat)),1)
plot(main="A) Mendel's ratios", bty='l', x=NULL, y=NULL,
xlab=expression(paste("Theoretical ", chi[1]^2)),
ylab=expression(paste("Mendel's ", chi[1]^2)),
xlim=c(0, 15), ylim=c(0, 15))
xint=seq(1/nrow(dat), 1-1/nrow(dat), length.out = nrow(dat))
for(i in 1:(nrow(dat)-1)) {
sd1=sqrt(i*(nrow(dat)-i+1)/((nrow(dat)+1)^2*(nrow(dat)+2)))
sd2=sqrt((i+1)*(nrow(dat)-(i+1)+1)/((nrow(dat)+1)^2*(nrow(dat)+2)))
polygon(x=qchisq(c(xint[i], xint[i], xint[i+1], xint[i+1]), 1),
y=qchisq(c(xint[i]-sd1, xint[i]+sd1, xint[i+1]+sd2,xint[i+1]-sd2), 1), col="grey90", border = "grey90")
}
abline(a=0, b=1, col="grey", lty=2)
points(x=xobs,
y=dat$MendelChi[od],
col=dat$col[od], cex=ifelse(dat$col[od] != "black", 1, 0.5),
pch=16)
dat$Eobs1=dat$Total*(dat$FisherRatio1/(dat$FisherRatio1+dat$FisherRatio2))
dat$Eobs2=dat$Total*(dat$FisherRatio2/(dat$FisherRatio1+dat$FisherRatio2))
dat$chi=(dat$observed1-dat$Eobs1)^2/dat$Eobs1+(dat$observed2-dat$Eobs2)^2/dat$Eobs2
od=order(dat$chi)
plot(main="B) Fisher's ratios", bty='l', x=NULL, y=NULL,
xlab=expression(paste("Theoretical ", chi[1]^2)),
ylab=expression(paste("Mendel's observed ", chi[1]^2)),
xlim=c(0, 15), ylim=c(0, 15))
xint=seq(1/nrow(dat), 1-1/nrow(dat), length.out = nrow(dat))
for(i in 1:(nrow(dat)-1)) {
sd1=sqrt(i*(nrow(dat)-i+1)/((nrow(dat)+1)^2*(nrow(dat)+2)))
sd2=sqrt((i+1)*(nrow(dat)-(i+1)+1)/((nrow(dat)+1)^2*(nrow(dat)+2)))
polygon(x=qchisq(c(xint[i], xint[i], xint[i+1], xint[i+1]), 1),
y=qchisq(c(xint[i]-sd1, xint[i]+sd1, xint[i+1]+sd2,xint[i+1]-sd2), 1), col="grey90", border = "grey90")
}
abline(a=0, b=1, col="grey", lty=2)
points(x=xobs,
y=dat$chi[od], col=dat$col[od], cex=ifelse(dat$col[od] != "black", 1, 0.5),
pch=16)
dat$SimuObs1=ceiling(rbinom(nrow(dat), dat$Total, dat$FisherRatio1/(dat$FisherRatio1+dat$FisherRatio2)))
dat$SimuObs2=dat$Total-dat$SimuObs1
dat$Schi=(dat$SimuObs1-dat$Eobs1)^2/dat$Eobs1+(dat$SimuObs2-dat$Eobs2)^2/dat$Eobs2
Sod=order(dat$Schi)
plot(main="C) Simulated ratios", bty='l', x=NULL, y=NULL,
xlab=expression(paste("Theoretical ", chi[1]^2)),
ylab=expression(paste("Simulation observed ", chi[1]^2)),
xlim=c(0, 15), ylim=c(0, 15))
xint=seq(1/nrow(dat), 1-1/nrow(dat), length.out = nrow(dat))
for(i in 1:(nrow(dat)-1)) {
sd1=sqrt(i*(nrow(dat)-i+1)/((nrow(dat)+1)^2*(nrow(dat)+2)))
sd2=sqrt((i+1)*(nrow(dat)-(i+1)+1)/((nrow(dat)+1)^2*(nrow(dat)+2)))
polygon(x=qchisq(c(xint[i], xint[i], xint[i+1], xint[i+1]), 1),
y=qchisq(c(xint[i]-sd1, xint[i]+sd1, xint[i+1]+sd2,xint[i+1]-sd2), 1), col="grey90", border = "grey90")
}
abline(a=0, b=1, col="grey", lty=2)
points(x=xobs,
y=dat$Schi[Sod], col=dat$col[Sod], cex=ifelse(dat$col[Sod] != "black", 1, 0.5),
pch=16)
plot(main="F3", dat$chi[dat$col!="black"], dat$MendelChi[dat$col!="black"], col=dat$col[dat$col!="black"], pch=16, bty="l",
xlab="Fisher's", ylab="Mendel's", xlim=c(0, 5), ylim=c(0, 5), cex=0.5)
abline(a=0, b=1, col="grey", lty=2)
