library(psych) data<-read.csv("ISC_total.csv") head(data) dataL<-data[data$Pop=="L",] dataZ<-data[data$Pop=="Z",] dataM<-data[data$Pop=="M",] dataS<-data[data$Pop=="S",] #L data pL<-c(dataL$P2,1-dataL$P2) pLgen1<-pL[1:16] pLgen2<-pL[17:32] pLgen3<-pL[33:48] pLgen4<-pL[49:64] pLgen5<-pL[65:96] pLgen6<-pL[97:128] op2L<-c(mean(pLgen1),mean(pLgen2),mean(pLgen3),mean(pLgen4)) op1L<-c(mean(pLgen5),mean(pLgen6)) vw<-(2/6)*var(op1L)+(4/6)*var(op2L) vb<-(mean(op1L)-mean(op2L))^2*(2/6)*(4/6) ovtL<-vw+vb #Zdata pZ<-c(dataZ$P2,1-dataZ$P2) pZgen1<-pZ[1:16] pZgen2<-pZ[17:32] pZgen3<-pZ[33:48] pZgen4<-pZ[49:64] pZgen5<-pZ[65:96] pZgen6<-pZ[97:128] op2Z<-c(mean(pZgen1),mean(pZgen2),mean(pZgen3),mean(pZgen4)) op1Z<-c(mean(pZgen5),mean(pZgen6)) vw<-(2/6)*var(op1Z)+(4/6)*var(op2Z) vb<-(mean(op1Z)-mean(op2Z))^2*(2/6)*(4/6) ovtZ<-vw+vb #Mdata pM<-c(dataM$P2,1-dataM$P2) pMgen1<-pM[1:16] pMgen2<-pM[17:32] pMgen3<-pM[33:48] pMgen4<-pM[49:64] pMgen5<-pM[65:96] pMgen6<-pM[97:128] op2M<-c(mean(pMgen1),mean(pMgen2),mean(pMgen3),mean(pMgen4)) op1M<-c(mean(pMgen5),mean(pMgen6)) vw<-(2/6)*var(op1M)+(4/6)*var(op2M) vb<-(mean(op1M)-mean(op2M))^2*(2/6)*(4/6) ovtM<-vw+vb #Sdata pS<-c(dataS$P2,1-dataS$P2) pSgen1<-pS[1:16] pSgen2<-pS[17:32] pSgen3<-pS[33:48] pSgen4<-pS[49:64] pSgen5<-pS[65:96] pSgen6<-pS[97:128] op2S<-c(mean(pSgen1),mean(pSgen2),mean(pSgen3),mean(pSgen4)) op1S<-c(mean(pSgen5),mean(pSgen6)) vw<-(2/6)*var(op1S)+(4/6)*var(op2S) vb<-(mean(op1S)-mean(op2S))^2*(2/6)*(4/6) ovtS<-vw+vb nsamp<-1000 varL<-rep() varZ<-rep() varM<-rep() varS<-rep() for (i in 1:nsamp){ #L loop p2_1<-mean(sample(pLgen1,16,replace=T)) p2_2<-mean(sample(pLgen2,16,replace=T)) p2_3<-mean(sample(pLgen3,16,replace=T)) p2_4<-mean(sample(pLgen4,16,replace=T)) p2<-c(p2_1,p2_2,p2_3,p2_4) p2L<-p2 p1_1<-mean(sample(pLgen5,16,replace=T)) p1_2<-mean(sample(pLgen6,16,replace=T)) p1<-c(p1_1,p1_2) vw<-(2/6)*var(p1)+(4/6)*var(p2) vb<-(mean(p1)-mean(p2))^2*(2/6)*(4/6) vtL<-vw+vb varL<-c(varL,vtL) #Z loop p2_1<-mean(sample(pZgen1,16,replace=T)) p2_2<-mean(sample(pZgen2,16,replace=T)) p2_3<-mean(sample(pZgen3,16,replace=T)) p2_4<-mean(sample(pZgen4,16,replace=T)) p2<-c(p2_1,p2_2,p2_3,p2_4) p2Z<-p2 p1_1<-mean(sample(pZgen5,16,replace=T)) p1_2<-mean(sample(pZgen6,16,replace=T)) p1<-c(p1_1,p1_2) vw<-(2/6)*var(p1)+(4/6)*var(p2) vb<-(mean(p1)-mean(p2))^2*(2/6)*(4/6) vtZ<-vw+vb varZ<-c(varZ,vtZ) # M loop p2_1<-mean(sample(pMgen1,16,replace=T)) p2_2<-mean(sample(pMgen2,16,replace=T)) p2_3<-mean(sample(pMgen3,16,replace=T)) p2_4<-mean(sample(pMgen4,16,replace=T)) p2<-c(p2_1,p2_2,p2_3,p2_4) p2M<-p2 p1_1<-mean(sample(pMgen5,16,replace=T)) p1_2<-mean(sample(pMgen6,16,replace=T)) p1<-c(p1_1,p1_2) vw<-(2/6)*var(p1)+(4/6)*var(p2) vb<-(mean(p1)-mean(p2))^2*(2/6)*(4/6) vtM<-vw+vb varM<-c(varM,vtM) # S loop p2_1<-mean(sample(pSgen1,16,replace=T)) p2_2<-mean(sample(pSgen2,16,replace=T)) p2_3<-mean(sample(pSgen3,16,replace=T)) p2_4<-mean(sample(pSgen4,16,replace=T)) p2<-c(p2_1,p2_2,p2_3,p2_4) p2S<-p2 p1_1<-mean(sample(pSgen5,16,replace=T)) p1_2<-mean(sample(pSgen6,16,replace=T)) p1<-c(p1_1,p1_2) vw<-(2/6)*var(p1)+(4/6)*var(p2) vb<-(mean(p1)-mean(p2))^2*(2/6)*(4/6) vtS<-vw+vb varS<-c(varS,vtS) } dL<-ovtL-varL dZ<-ovtZ-varZ dM<-ovtM-varM dS<-ovtS-varS #L confidence interval ovtL-quantile(dL,0.975) ovtL-quantile(dL,0.025) #Z confidence interval ovtZ-quantile(dZ,0.975) ovtZ-quantile(dZ,0.025) #M confidence interval ovtM-quantile(dM,0.975) ovtM-quantile(dM,0.025) #S confidence interval ovtS-quantile(dS,0.975) ovtS-quantile(dS,0.025)