pca_log <- prcomp(log(iris[,1:4]),center = TRUE,scale. = TRUE) # print(pca_log) install.packages("ggbiplot") ?prcomp pca1<-prcomp(pca_log,scale.=TRUE) plot(pca_log) plot(pca_log$rotation[,1],pca_log$rotation[,2],col="red") Allozyme_data<-matrix(NA, nrow=9,ncol=20) # # Lap1<-c(0,0,0,0,0,0,0,0.025,0.030)# Lap2<-c(0,0,0,0,0,0,0,0,0.030)# Lap3<-c(0,0,0,0,0,0,0,0.350,0)# Lap4<-c(0,0.214,0.250,0,0.140,0,0.024,0.475,0.156)# Lap5<-c(0.167,0.119,0.286,0.278,0.160,0.136,0.119,0,0.345)# Lap6<-c(0.833,0.667,0.214,0.611,0.580,0.705,0.666,0,0.345)# Lap7<-c(0,0,0.250,0.111,0.120,0.091,0.167,0.150,0.094)# Lap8<-c(0,0,0,0,0,0.068,0.024,0,0)# # Mpi1<-c(0,0.024,0,0,0.040,0.023,0.050,0,0)# Mpi2<-c(0.167,0,0.107,0,0.060,0.068,0.025,0,0)# Mpi3<-c(0.833,0.262,0.107,0.056,0.360,0.045,0.125,0.263,0.500)# Mpi4<-c(0,0.547,0.607,0.721,0.460,0.546,0.650,0.737,0.406)# Mpi5<-c(0,0.143,0.143,0.167,0.080,0.273,0.150,0,0.094)# Mpi6<-c(0,0.024,0.036,0.056,0,0.045,0,0,0)# # Pgi1<-c(0,0,0.036,0,0.040,0,0.190,0.265,0.056)# Pgi2<-c(0,0,0,0,0,0.023,0,0,0)# Pgi3<-c(0.667,0.095,0.214,0.111,0.140,0.931,0.525,0.059,0.167)# Pgi4<-c(0.333,0.429,0.214,0.278,0.240,0.023,0.190,0.529,0.527)# Pgi5<-c(0,0.309,0.429,0.444,0.480,0.023,0.095,0.147,0.250)# Pgi6<-c(0,0.167,0.107,0.167,0.100,0,0,0,0)# Allozyme_data[,1]<-Lap1# Allozyme_data[,2]<-Lap2# Allozyme_data[,3]<-Lap3# Allozyme_data[,4]<-Lap4# Allozyme_data[,5]<-Lap5# Allozyme_data[,6]<-Lap6# Allozyme_data[,7]<-Lap7# Allozyme_data[,8]<-Lap8# Allozyme_data[,9]<-Mpi1# Allozyme_data[,10]<-Mpi2# Allozyme_data[,11]<-Mpi3# Allozyme_data[,12]<-Mpi4# Allozyme_data[,13]<-Mpi5# Allozyme_data[,14]<-Mpi6# Allozyme_data[,15]<-Pgi1# Allozyme_data[,16]<-Pgi2# Allozyme_data[,17]<-Pgi3# Allozyme_data[,18]<-Pgi4# Allozyme_data[,19]<-Pgi5# Allozyme_data[,20]<-Pgi6# # rownames <-c("Trat", "Chon_Buri", "Chumphon_5", "Chumphon_6", "Surat_Thani", "Satun", "Trang", "Sydney_Aust", "Magnetic_Island")# colnames(Allozyme_data)<-c( "Lap1", "Lap2", "Lap3", "Lap4","Lap5", "Lap6"," Lap7","Lap8", "Mpi1","Mpi2","Mpi3", "Mpi4", "Mpi5","Mpi6", "Pgi1","Pgi2","Pgi3","Pgi4","Pgi5","Pgi6")# # pca1<-prcomp(t(Allozyme_data), scale. = TRUE)# head(pca1$rotation,3)# colours1<-c("orange","orange","orange","orange","orange","orange", "orange", "blue","blue")# # par(mfrow=c(2,2))# plot(pca1$rotation[,1], pca1$rotation[,2], col=colours1, pch=15)# plot(pca1$rotation[,2], pca1$rotation[,3], col=colours1, pch=15)# plot(pca1$rotation[,1], pca1$rotation[,3], col=colours1, pch=15)# plot(pca1$rotation[,2], pca1$rotation[,4], col=colours1, pch=15) par(mfrow=c(1,1))# localities<-c("Trat", "Chon_Buri", "Chumphon_5", "Chumphon_6", "Surat_Thani", "Satun", "Trang", "Sydney_Aust", "Magnetic_Island")# plot(pca1$rotation[,1], pca1$rotation[,2], col=colours1, pch=15)# xoffset<-0.01*(max(pca1$rotation[,1])-min(pca1$rotation[,1]))# yoffset<-0.02*(max(pca1$rotation[,2])-min(pca1$rotation[,2]))# text(pca1$rotation[,1]+xoffset, pca1$rotation[,2]+yoffset,localities,cex=0.7,xpd=NA) head(pca1) rownames(pca_log)<-iris[,5] rownames(newiris)<-iris[,5] newiris<-iris# rownames(newiris)<-iris[,5] class(newiris[,5]) newiris[,5]<-as.numeric(newiris[,5] ) newiris lognewiris[,1:4]<-log(newiris[,1:4]) lognewiris<-newiris# lognewiris[,1:4]<-log(newiris[,1:4]) lognewiris pca_log <- prcomp(lognewiris,center = TRUE,scale. = TRUE) print(pca_log) summary(pca_log) library(devtools)# install_github("ggbiplot", "vqv")# library(ggbiplot) install_github("vqv/ggbiplot", force = TRUE) print(gglogirispcs) ?ggbiplot cov(pca_log[,1…4]) cov(pca_log[,1:4]) cov(iris[,1:4]) library(ggfortify) install.packages("ggfortify") data(iris) # iris is in Base R# newiris<-iris PCA1<-prcomp(newiris[,1:4],center = TRUE,scale. = TRUE) pca2_plot<-autoplot(PCA1,data=newiris,colour='Species')# pca2_plot ?ggfortify ??ggfortify library(ggplot2) ?ggfortify pca2_plot<-autoplot(PCA1,data=newiris,colour='Species')# pca2_plot library(ggfortify) pca2_plot<-autoplot(PCA1,data=newiris,colour='Species')# pca2_plot ?pccomp ?prcomp prcomp(USArrests, scale = TRUE)# prcomp(~ Murder + Assault + Rape, data = USArrests, scale = TRUE)# plot(prcomp(USArrests))# summary(prcomp(USArrests, scale = TRUE))# biplot(prcomp(USArrests, scale = TRUE)) par(mfrow=c(2,2))# #### basic PCA on raw data# PCA1<-prcomp(newiris[,1:4],center = TRUE,scale. = TRUE) # plot(PCA1)# #### with logged data# PCA_log <- prcomp(log(newiris[,1:4]),center = TRUE,scale. = TRUE) # plot(PCA_log) biplot(prcomp(PCA1, scale = TRUE)) biplot(PCA1, scale = TRUE)) biplot(PCA1, scale = TRUE) ?biplot biplot(pca1) biplot(pca1,col=c(1:9)) ?biplot pca1 biplot(pca1[,2:3],col=c(1:9)) biplot(pca1[,2],pca1[,3],col=c(1:9)) biplot(princomp(USArrests)) biplot(pca1) biplot(princomp(Allozyme_data)) biplot.princomp(Allozyme_data) biplot(prcomp(USArrests, scale = TRUE)) plot(pca1[,2],paca1[,3]) pca1 plot(pca1$PC2,paca1$PC3) plot(pca1$PC2,pca1$PC3) pca1$PC3 summary(pca1) biplot(pca1) par(mfrow=c(2,2)) # PCA1<-prcomp(iris[,1:4],center = TRUE,scale. = TRUE) # plot(PCA1) # PCA_log <- prcomp(log(iris[,1:4]),center = TRUE,scale. = TRUE) # plot(PCA_log) PCA1<-prcomp(iris[,1:4],center = TRUE,scale. = TRUE) # biplot(PCA1) # PCA_log <- prcomp(log(iris[,1:4]),center = TRUE,scale. = TRUE) # biplot(PCA_log) PCA1<-prcomp(iris[,1:4],center = TRUE,scale. =TRUE) # biplot(PCA1,col=as.numeric(iris[,5])) # PCA_log <- prcomp(log(iris[,1:4]),center = TRUE,scale. =TRUE) # biplot(PCA_log, col=as.numeric(iris[,5]))) ?biplot PCA1<-prcomp(iris[,1:4],center = TRUE,scale. =TRUE) # biplot(PCA1,col=as.numeric(iris[,5])) # PCA_log <- prcomp(log(iris[,1:4]),center = TRUE,scale. =TRUE) # biplot(PCA_log, col=c("red,"blue)) PCA1<-prcomp(iris[,1:4],center = TRUE,scale. =TRUE) # biplot(PCA1,col=as.numeric(iris[,5])) # PCA_log <- prcomp(log(iris[,1:4]),center = TRUE,scale. =TRUE) # biplot(PCA_log, col=c("red","blue)) ) biplot(PCA_log,col=c("red","blue)) biplot(PCA_log,col=c("red","blue")) biplot(PCA_log,col=c("red","blue"),cex=0.5) par(mfrow=c(2,2)) # PCA1<-prcomp(iris[,1:4],center = TRUE,scale. =TRUE) # biplot(PCA1) # PCA_log <- prcomp(log(iris[,1:4]),center = TRUE,scale. =TRUE) # biplot(PCA_log) summary(PCA1) summary(PCA_LOG) summary(PCA_log) yv<-predict(PCA_log[,1]) yv<-predict(PCA_log)[,1] yv2<-predict(PCA_log)[,2] pca2_plot<-autoplot(PCA1,data=iris,colour='Species') # pca2_plot ?autoplot pca2_plot<-autoplot(PCA1,data=iris,colour='Species') # pca2_plot pca2_plot<-autoplot(PCA1,data=iris,colour='Species') pca2_plot ?slplot ??slplot data.pca$pca_log data.pca_log$pca_log summary(PCA_log) PCA_log$PC1 head(iris) plot(PCA1$iris[,1],PCA1$iris[,2]) PCA1 plot(PCA1$PC1,PCA1$PC2) PCA1$PC1 PCA1$iris PCA1$iris[,1] atr(PCA1) attr(PCA1) str(PCA1) loadings <- PCA1$rotation loadings axes <- predict(PCA1, newdata =Sepal.Length) pca <- prcomp(~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, data = iris, scale = TRUE) pca (loadings <- pca$rotation) axes <- predict(pca, newdata = USArrests)# head(axes, 4) axes <- predict(pca, newdata = iris) axes axes <- predict(PCA1, newdata = iris) irisPCAs <- predict(PCA1, newdata = iris) tail(iris) class(irisPCAs) head(irisPCAs) colours1<-5+as.numeric(iris[,5])# par(mfrow=c(2,2))# plot(irisPCAs$PC1,irisPCAs$PC2,col=colours1) irisPCAs<-as.data.frame(irisPCAs) colours1<-5+as.numeric(iris[,5])# par(mfrow=c(2,2))# plot(irisPCAs$PC1,irisPCAs$PC2,col=colours1) plot(irisPCAs$PC1,irisPCAs$PC3,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC4,col=colours1)# plot(irisPCAs$PC2,irisPCAs$PC3,col=colours1) colours1<-as.numeric(iris[,5])+8# par(mfrow=c(2,2))# plot(irisPCAs$PC1,irisPCAs$PC2,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC3,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC4,col=colours1)# plot(irisPCAs$PC2,irisPCAs$PC3,col=colours1) colours1<-as.numeric(iris[,5])+9# par(mfrow=c(2,2))# plot(irisPCAs$PC1,irisPCAs$PC2,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC3,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC4,col=colours1)# plot(irisPCAs$PC2,irisPCAs$PC3,col=colours1) levels(iris[,5]) ?legensd ?legend legend(-6,2,xpd=NA,levels(iris[,5]) pch=1, col=as.numeric(levels(iris[,5])+9)+ legend(-6,2,xpd=NA,levels(iris[,5]) pch=1, col=as.numeric(levels(iris[,5])+9) legend(-6,2,xpd=NA,levels(iris[,5]), pch=1, col=as.numeric(levels(iris[,5])+9) ) as.numeric(levels(iris[,5]) ) levels(iris[,5]) as.numeric(levels(iris[,5])) legend(-6,2,xpd=NA,c("setosa","versicolor","virginica"), pch=1, col=c("red","green","blue")) par(mfrow=c(2,2))# plot(irisPCAs$PC1,irisPCAs$PC2,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC3,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC4,col=colours1)# plot(irisPCAs$PC2,irisPCAs$PC3,col=colours1)# # legend(-5,2,xpd=NA,c("setosa","versicolor","virginica"), pch=1, col=c("red","green","blue")) par(mfrow=c(2,2))# plot(irisPCAs$PC1,irisPCAs$PC2,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC3,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC4,col=colours1)# plot(irisPCAs$PC2,irisPCAs$PC3,col=colours1)# # legend(-5.2,2,xpd=NA,c("setosa","versicolor","virginica"), pch=1, col=c("red","green","blue")) par(mfrow=c(2,2))# plot(irisPCAs$PC1,irisPCAs$PC2,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC3,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC4,col=colours1)# plot(irisPCAs$PC2,irisPCAs$PC3,col=colours1)# # legend(-5.15,2,xpd=NA,c("setosa","versicolor","virginica"), pch=1, col=c("red","green","blue")) par(mfrow=c(2,2))# plot(irisPCAs$PC1,irisPCAs$PC2,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC3,col=colours1)# plot(irisPCAs$PC1,irisPCAs$PC4,col=colours1)# plot(irisPCAs$PC2,irisPCAs$PC3,col=colours1)# # legend(-5.15,1.9,xpd=NA,c("setosa","versicolor","virginica"), pch=1, col=c("red","green","blue")) library(rworldmap)# data(countryExData)# country2Region(regionType="Stern") # sternEnvHealth<-country2Region(inFile=countryExData, nameDataColumn="ENVHEALTH",joinCode="ISO3", nameJoinColumn="ISO3V10",regionType="Stern",FUN="mean")# sPDF<-joinCountryData2Map(countryExData,joinCode = "ISO3",nameJoinColumn="ISO3V10")# mapCountryData(sPDF, nameColumnToPlot="BIODIVERSITY", mapTitle="World Biodiversity", oceanCol="lightblue", missingCountryCol="white", borderCol="black")# ####### mapCountryData(sPDF,nameColumnToPlot="BIODIVERSITY", mapTitle="Asia",oceanCol="lightblue",missingCountryCol="white", mapRegion="Asia",borderCol="black")# ####### mapCountryData(sPDF,nameColumnToPlot="BIODIVERSITY", mapTitle="Australia",oceanCol="lightblue",missingCountryCol="white", xlim=c(110,160), ylim=c(-48,-5),borderCol="black")# # mapBubbles(dF=getMap())# mapBubbles(dF=getMap(),oceanCol="lightblue",landCol="wheat")# head(getMap()$NAME)# # mammals<- read.table(text="# Country 'Mammal species' 'threatened spp' Lat Lon# Canada 426 40 60 -95# Brazil 648 80 -34 -64# Mexico 923 96 24 -102# Peru 647 63 -9 -75# Colombia 442 58 3.9 -73# 'United States' 440 40 39 -99# Argentina 374 38 -35 -65# Ecuador 372 47 -1 -78# Bolivia 362 21 -17 -64# Chile 143 19 -37 -72# ", header=TRUE) # # mapBubbles(dF=mammals, nameZSize="Mammal.species",nameZColour="Country",oceanCol="lightblue", landCol="wheat",nameX = "Lat", nameY = "Lon", main="Number of Mammal species",symbolSize=0.5) data(countryExData)# country2Region(regionType="Stern") # sternEnvHealth<-country2Region(inFile=countryExData, nameDataColumn="ENVHEALTH",joinCode="ISO3", nameJoinColumn="ISO3V10",regionType="Stern",FUN="mean")# sPDF<-joinCountryData2Map(countryExData,joinCode = "ISO3",nameJoinColumn="ISO3V10")# mapCountryData(sPDF, nameColumnToPlot="BIODIVERSITY", mapTitle="World Biodiversity", oceanCol="lightblue", missingCountryCol="white", borderCol="black") mosses<-c(53, 85, 69, 98, 130, 90)# liverworts<-c(47, 67, 85, 113, 190, 239)# data<-rbind(mosses,liverworts)# barplot(data,beside=TRUE) x<-rchisq(sample_size,5) # hist(x,prob=T,col="grey30",ylab="Probability") # curve(dchisq(x,df=2),lwd=2,col="green",add=TRUE) # curve(dchisq(x,df=5),lwd=2,col="blue",add=TRUE) # curve(dchisq(x,df=10),lwd=2,col="red",add=TRUE) # leg_text<-c("d.f. = 2","d.f. = 5","d.f.=10")# legend("topright",lwd=2,col=c("green","blue","red"),leg_text) hist(x,prob=T,col="grey30",ylab="Probability") # curve(dchisq(x,df=2),lwd=2,col="green",add=TRUE) # curve(dchisq(x,df=5),lwd=2,col="blue",add=TRUE) # curve(dchisq(x,df=10),lwd=2,col="red",add=TRUE) # leg_text<-c("d.f. = 2","d.f. = 5","d.f.=10") # legend("topright",lwd=2,col=c("green","blue","red"),leg_text) sample_size<-10000 # x<-rchisq(sample_size,5) # hist(x,col="grey30",ylab="Frequency") curve(dchisq(x,df=2),lwd=2,col="green",add=TRUE) # curve(dchisq(x,df=5),lwd=2,col="blue",add=TRUE) # curve(dchisq(x,df=10),lwd=2,col="red",add=TRUE) # leg_text<-c("d.f. = 2","d.f. = 5","d.f.=10") # legend("topright",lwd=2,col=c("green","blue","red"),leg_text) sample_size<-10000# x<-rchisq(sample_size,5) # hist(x,col="grey30",ylab="Frequency") # curve(dchisq(x,df=2),lwd=2,col="green",add=TRUE) # curve(dchisq(x,df=5),lwd=2,col="blue",add=TRUE) # curve(dchisq(x,df=10),lwd=2,col="red",add=TRUE) # leg_text<-c("d.f. = 2","d.f. = 5","d.f.=10")# legend("topright",lwd=2,col=c("green","blue","red"),leg_text) sample_size<-10000# x<-rchisq(sample_size,5) # hist(x,prob=T,col="grey30",ylab="Probability") # curve(dchisq(x,df=2),lwd=2,col="green",add=TRUE) # curve(dchisq(x,df=5),lwd=2,col="blue",add=TRUE) # curve(dchisq(x,df=10),lwd=2,col="red",add=TRUE) # leg_text<-c("d.f. = 2","d.f. = 5","d.f.=10")# legend("topright",lwd=2,col=c("green","blue","red"),leg_text) sample_size<-10000# x<-rchisq(sample_size,5) # hist(x,prob=TRUE,col="grey30",ylab="Probability") # curve(dchisq(x,df=2),lwd=2,col="green",add=TRUE) # curve(dchisq(x,df=5),lwd=2,col="blue",add=TRUE) # curve(dchisq(x,df=10),lwd=2,col="red",add=TRUE) # leg_text<-c("d.f. = 2","d.f. = 5","d.f.=10")# legend("topright",lwd=2,col=c("green","blue","red"),leg_text) mosses<-c(53, 85, 69, 98, 130, 90)# liverworts<-c(47, 67, 85, 113, 190, 239)# data<-rbind(mosses,liverworts)# # altitude<-c(400,600,800,1000,1200, 1300)# colnames(data)<- altitude# barplot(data, col=c("lightgreen","darkgreen"),beside=TRUE, xlab="Altitude (m)",ylab="Number of bryophytes")# legend("topleft", fill=c("lightgreen","darkgreen"), c("mosses","liverworts"), cex=1) #Example 3 - Spatial host parasitoid model# # Z<-30# H<- matrix (NA, nrow=Z, ncol=Z) # P<- matrix (NA, nrow=Z, ncol=Z) # # G<- 100 # # Beddington<-function(fNt,fPt,fr,falpha,fK,fc){# NtPlus1<-fNt*exp((fr*(1-(fNt/fK)))-(falpha*fPt))# PtPlus1<-fc*fNt*(1-(exp(-falpha*fPt)))# answer<-c(NtPlus1,PtPlus1)# answer# }# # random_asign_to_neighbours<-function(to_Add_f, Number_moving, neighbours){# cell<-rep(NA,8)# breaks<- sort(sample(1: Number_moving, 7, replace=TRUE))# for(i in 2:7) {cell[i] <-breaks[i]- breaks[i-1]}# cell[1]<-breaks[1]# cell[8]<- Number_moving - breaks[7]# for(x in 1:8){# to_Add_f[neighbours[x,1],neighbours [x,2]]<- to_Add_f[neighbours[x,1],neighbours [x,2]] + cell[x]# } # to_Add_f# } # random_asign_to_neighbours function# H[1:Z,1:Z]<-100# P[1:Z,1:Z]<-10# alpha<-0.0025# r<-2.1# c<-1# K<-1000 # proportion_to_migrate<-0.1# neighbours<-matrix(NA, nrow=8, ncol=2)# par(mfrow=c(3,2))# save_runs<-c(10,30,60,100) # hold_H<-matrix(NA, nrow=Z*Z, ncol=length(save_runs))# 4 runs# hold_P<- matrix(NA, nrow=Z*Z, ncol=length(save_runs))# count<-0# # for(g in 1:G) { # for(ba in 1:Z){ # for(bb in 1:Z){# answer<-Beddington(H[ba,bb],P[ba,bb],r,alpha,K,c)# H[ba,bb]<-round(answer[1])# P[ba,bb]<-round(answer[2])# } # for bb loop# } # for ba loop# tomove_Host<-round(H* proportion_to_migrate)# tomove_Parasitoid<-round(P* proportion_to_migrate)# H<- round(H-tomove_Host) # P<- round(P-tomove_Parasitoid) # to_add_Host<-matrix(0, nrow=Z, ncol=Z)# to_add_Parasitoid<-matrix(0, nrow=Z, ncol=Z)# # for(ba in 1:Z){# for(bb in 1:Z){ # neighbours[1,1:2]<- c(ba-1, bb-1) # neighbours[2,1:2]<- c(ba-1, bb)# neighbours[3,1:2]<- c(ba-1, bb+1)# neighbours[4,1:2]<- c(ba, bb-1)# neighbours[5,1:2]<- c(ba, bb+1)# neighbours[6,1:2]<- c(ba+1,bb-1)# neighbours[7,1:2]<- c(ba+1,bb)# neighbours[8,1:2]<- c(ba+1,bb+1)# neighbours[which(neighbours==0)]<- Z # neighbours[which(neighbours==Z+1)]<-1 # to_add_Host <- random_asign_to_neighbours (to_add_Host, H[ba,bb], neighbours)# to_add_Parasitoid <- random_asign_to_neighbours (to_add_Parasitoid, P[ba,bb], neighbours) # } #end bb loop# } # end ba loop# H<-H+ to_add_Host# P<- P+ to_add_Parasitoid# # if(g %in% save_runs){# count<-count+1# hold_H[,count]<- unlist(expand.grid(H))# hold_P[,count]<- unlist(expand.grid(P))# } #end if g in save list# } # for g loop# plot_expandgrid<-function(X,Z){ # mcol<-max(X)# start<-0# for(j in 1:Z){# for(k in 1:Z){# start<-start+1# points(j,k, col=round(10*X[start]/mcol),pch=15)# } # for k loop# } # for j loop# } # end plot_generation function# getcoloursrainbow<-function(x){# purplef<-220# f<-unlist(as.list(rainbow(256))) # mn<-0# mx<-max(x,na.rm=TRUE)# cor<-(x-mn)/(mx-mn)# ans<-f[(1+ceiling(purplef*(cor)))]# ans# }# par(mfrow=c(2,2))# # point_colours<- getcoloursrainbow(hold_H[,1])# point_colours2<-matrix(point_colours,nrow=Z)# plot(NULL, xlim=c(1,Z),ylim=c(1,Z), axes=F, xlab="",ylab="", main="Host")# for(j in 1:Z){# for(k in 1:Z){# points(j,k, col= point_colours2[j,k],pch=15)# }}# point_colours3<- getcoloursrainbow(hold_P[,1])# point_colours4<-matrix(point_colours3,nrow=Z)# plot(NULL, xlim=c(1,Z), ylim=c(1,Z), axes=F, xlab="" ,ylab="", main="Parasitoid")# for(j in 1:Z){# for(k in 1:Z){# points(j,k, col= point_colours4[j,k],pch=15)# }}# point_colours<- getcoloursrainbow(hold_H[,4])# point_colours2<-matrix(point_colours,nrow=Z)# plot(NULL, xlim=c(1,Z),ylim=c(1,Z), axes=F, xlab="",ylab="", main="Host")# for(j in 1:Z){# for(k in 1:Z){# points(j,k, col= point_colours2[j,k],pch=15)# }}# point_colours3<- getcoloursrainbow(hold_P[,4])# point_colours4<-matrix(point_colours3,nrow=Z)# plot(NULL, xlim=c(1,Z), ylim=c(1,Z), axes=F, xlab="" ,ylab="", main="Parasitoid")# for(j in 1:Z){# for(k in 1:Z){# points(j,k, col= point_colours4[j,k],pch=15)# }}# centreX<- -7# centreY<- 45# LX<-8# lowY<- centreY-2# highY<- centreY+2# increment<-0.08# line_colour<- getcoloursrainbow(1:220)# for(i in 1:220){# X<- centreX-LX+(i*increment)# lines (c(X,X),c(lowY,highY), col=line_colour[i], xpd=NA)# }# text(xpd=NA, -24, centreY, "Population = 0")# text(xpd=NA, 15, centreY, "Maximum population")# # text(centreX, centreY+6,paste("After",save_runs[1],"generations", sep=" "),xpd=NA)# text(centreX, -5,paste("After",save_runs[4],"generations", sep=" "),xpd=NA) library(MASS)# data(menarche)# with(menarche,plot(Age, Menarche/Total,pch=19,col="blue"))# Mprob <- glm(cbind(Menarche,Total - Menarche) ~ Age,binomial(link = probit),data = menarche)# xv<-seq(min(menarche$Age),max(menarche$Age),length.out=500)# yv<-predict(Mprob,list(Age =xv),type="response")# lines(xv,yv,lwd=2,col="red")# dose.p(Mprob,cf=c(1,2),p=c(0.1,0.5,0.9)) HW<-function(a1,a2,a3,a4,a5){# x<-list(NULL)# x<-c(x,a1^2)# x<-c(x,2*a1*a2)# x<-c(x,2*a1*a3)# x<-c(x,2*a1*a4)# x<-c(x,2*a1*a5)# x<-c(x,a2^2)# x<-c(x,2*a2*a3)# x<-c(x,2*a2*a4)# x<-c(x,2*a2*a5)# x<-c(x,a3^2)# x<-c(x,a3*a4)# x<-c(x,a3*a5)# x<-c(x,a4^2)# x<-c(x,a4*a5)# x<-c(x,a5^2)# return(x)# }# # start_allele<-c(0.2,0.2,0.2,0.2,0.2)# # HW(start_allele) HW(start_allele[1:5]) HW(start_allele[1],start_allele[2],start_allele[3],start_allele[4],start_allele[5]) HW<-function(a1,a2,a3,a4,a5){# x<-list(NULL)# x<-c(x,a1^2)# x<-c(x,2*a1*a2)# x<-c(x,2*a1*a3)# x<-c(x,2*a1*a4)# x<-c(x,2*a1*a5)# x<-c(x,a2^2)# x<-c(x,2*a2*a3)# x<-c(x,2*a2*a4)# x<-c(x,2*a2*a5)# x<-c(x,a3^2)# x<-c(x,a3*a4)# x<-c(x,a3*a5)# x<-c(x,a4^2)# x<-c(x,a4*a5)# x<-c(x,a5^2)# return(unlist(x))# } HW(start_allele[1],start_allele[2],start_allele[3],start_allele[4],start_allele[5]) hw<-HW(start_allele[1],start_allele[2],start_allele[3],start_allele[4],start_allele[5]) sum(hw) HW<-function(a1,a2,a3,a4,a5){# x<-list(NULL)# x<-c(x,a1^2)# x<-c(x,2*a1*a2)# x<-c(x,2*a1*a3)# x<-c(x,2*a1*a4)# x<-c(x,2*a1*a5)# x<-c(x,a2^2)# x<-c(x,2*a2*a3)# x<-c(x,2*a2*a4)# x<-c(x,2*a2*a5)# x<-c(x,a3^2)# x<-c(x,2*a3*a4)# x<-c(x,2*a3*a5)# x<-c(x,a4^2)# x<-c(x,2*a4*a5)# x<-c(x,a5^2)# return(unlist(x))# }# # start_allele<-c(0.2,0.2,0.2,0.2,0.2)# # hw<-HW(start_allele[1],start_allele[2],start_allele[3],start_allele[4],start_allele[5]) sum(hw) hw<-1000*HW(start_allele[1],start_allele[2],start_allele[3],start_allele[4],start_allele[5]) hw sum(hw) knor<-read.table("krill.txt",header=TRUE) par(mfrow=c(2,2)) # plot(knor$wwt_grams,knor$ulO2_min,ylab="Oxygen consumption (ul/min)",xlab="Wet weight(grams)") # plot(log10(knor$wwt_grams),log10(knor$ulO2_min),ylab="log(Oxygen consumption (ul/min))",xlab="log(Wet weight(grams))") par(mfrow=c(2,2)) # model1<-lm(log(knor$ulO2_min) ~log(knor$wwt_grams)) # plot(model1) model3I<-lm(log(knor$ulO2_min) ~ knor$wwt_grams +I(knor$wwt_grams ^2)+ I(knor$wwt_grams^3)) model2<-lm(log(knor$ulO2_min) ~ poly(log(knor$wwt_grams),2)) # anova(model2,model3b) anova(model2,modelpoly3) modelpoly3<-lm(log(knor$ulO2_min) ~ poly(log(knor$wwt_grams),3)) anova(model2,modelpoly3) cerana<-c(1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,0,1,0,1,1,1,0,1,1,0,1,1,0,0,0,1,0,1,0,0,0,0,1, 1,0,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,1,0,0,1,1,0,0,0,1,0,1, 1,0,0,1,1,1,0,1,1,1,1,0,1,1,1) # dorsata<-c(0,1,1,1,0,0,1,1,0,0,0,0,0,1,0,1,1,1,0,0,1,0,0,0,0,1,0,1,0,0,1,0,0,1,1,1,1,0, 1,0,1,1,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,1,0,0,0,0,1,0, 0,0,0,1,0,0,0,1,0,1,0,1,0,0,1,0) # florea<-c(0,0,0,1,1,1,0,1,0,1,1,1,0,0,1,0,0,0,0,1,1,0,1,1,0,0,0,1,1,1,0,0,0,1,0,0,0,0,1, 0,0,1,1,0,1,1,0,0,1,0,0,0,1,1,0,0,0,1,1,1,1,1,1,1,0,0,0,1,1,0,0,0,0,1,1,0,1,0,0,1,0,1,0,0,1,1,0,0,1,1,1,0,0,1,0) # mellifera<-c(0,1,0,1,0,0,0,1,1,1,0,1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,1,0,0,0,0,1,1,1,1,0,0,0, 1,1,0,0,1,0,1,1,1,0,1,0,0,1,1,1,0,0,1,0,0,1,0,1,1,1,1,1,1,0,0,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,0,0,0,0,1,1) # data<-cbind(cerana,dorsata,florea,mellifera) cerana<-c(1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,0,1,0,1,1,1,0,1,1,0,1,1,0,0,0,1,0,1,0,0,0,0,1,1,0,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,1,0,0,1,1,0,0,0,1,0,1,1,0,0,1,1,1,0,1,1,1,1,0,1,1,1)# dorsata<-c(0,1,1,1,0,0,1,1,0,0,0,0,0,1,0,1,1,1,0,0,1,0,0,0,0,1,0,1,0,0,1,0, 0,1,1,1,1,0,1,0,1,1,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,1,0,1,0,1,0,0,1,0)# florea<-c(0,0,0,1,1,1,0,1,0,1,1,1,0,0,1,0,0,0,0,1,1,0,1,1,0,0,0,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,1,0,1,1,0,0,1,0,0,0,1,1,0,0,0,1,1,1,1,1,1,1,0,0,0,1,1,0,0,0,0,1,1,0,1,0,0,1,0,1,0,0,1,1,0,0,1,1,1,0,0,1,0)# mellifera<-c(0,1,0,1,0,0,0,1,1,1,0,1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,1,0,0,0,0,1,1,1,1,0,0,01,1,0,0,1,0,1,1,1,0,1,0,0,1,1,1,0,0,1,0,0,1,0,1,1,1,1,1,1,0,0,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,0,0,0,0,1,1)# data<-cbind(cerana,dorsata,florea,mellifera) cerana<-c(1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,0,1,0,1,1,1,0,1,1,0,1,1,0,0,0,1,0,1,0,0,0,0,1, 1,0,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,1,0,0,1,1,0,0,0,1,0,1, 1,0,0,1,1,1,0,1,1,1,1,0,1,1,1) # dorsata<-c(0,1,1,1,0,0,1,1,0,0,0,0,0,1,0,1,1,1,0,0,1,0,0,0,0,1,0,1,0,0,1,0, 0,1,1,1,1,0, 1,0,1,1,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,1,0,0,0,0,1,0, 0,0,0,1,0,0,0,1,0,1,0,1,0,0,1,0) # florea<-c(0,0,0,1,1,1,0,1,0,1,1,1,0,0,1,0,0,0,0,1,1,0,1,1,0,0,0,1,1,1,0,0,0,1,0,0,0,0,1, 0,0,1,1,0,1,1,0,0,1,0,0,0,1,1,0,0,0,1,1,1,1,1,1,1,0,0,0,1,1,0,0,0,0,1,1,0,1,0,0,1, 0,1,0,0,1,1,0,0,1,1,1,0,0,1,0) # mellifera<-c(0,1,0,1,0,0,0,1,1,1,0,1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,1,0,0,0,0,1,1,1,1,0,0, 0,1,1,0,0,1,0,1,1,1,0,1,0,0,1,1,1,0,0,1,0,0,1,0,1,1,1,1,1,1,0,0,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,0,0,0,0,1,1) # data<-cbind(cerana,dorsata,florea,mellifera) lengths(cerana,dorsata,florea,mellifera) length(cerana,dorsata,florea,mellifera) length(cerana) length(dorsata) length(florea) length(mellifera) mellifera<-c(0,1,0,1,0,0,0,1,1,1,0,1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,1,0,0,0,0,1,1,1,1,0,0,01,1,0,0,1,0,1,1,1,0,1,0,0,1,1,1,0,0,1,0,0,1,0,1,1,1,1,1,1,0,0,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,0,0,0,0,1,1) length(mellifera) cerana<-c(1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,0,1,0,1,1,1,0,1,1,0,1,1,0,0,0,1,0,1,0,0,0,0,1,1,0,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,1,0,0,1,1,0,0,0,1,0,1,1,0,0,1,1,1,0,1,1,1,1,0,1,1,1)# dorsata<-c(0,1,1,1,0,0,1,1,0,0,0,0,0,1,0,1,1,1,0,0,1,0,0,0,0,1,0,1,0,0,1,0, 0,1,1,1,1,0,1,0,1,1,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,1,0,1,0,1,0,0,1,0)# florea<-c(0,0,0,1,1,1,0,1,0,1,1,1,0,0,1,0,0,0,0,1,1,0,1,1,0,0,0,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,1,0,1,1,0,0,1,0,0,0,1,1,0,0,0,1,1,1,1,1,1,1,0,0,0,1,1,0,0,0,0,1,1,0,1,0,0,1,0,1,0,0,1,1,0,0,1,1,1,0,0,1,0)# mellifera<-c(0,1,0,1,0,0,0,1,1,1,0,1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,1,0,0,0,0,1,1,1,1,0,0,1,1,0,0,1,0,1,1,1,0,1,0,0,1,1,1,0,0,1,0,0,1,0,1,1,1,1,1,1,0,0,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,0,0,0,0,1,1)# data<-cbind(cerana,dorsata,florea,mellifera)# raw_score<-sum(rowSums(data[,1:4])^2)# raw_score# n<-nrow(data)# randomisations<-NULL# for(i in 1:100000){# d1<-sample(data[,1],n,replace=F)# d2<-sample(data[,2],n,replace=F)# d3<-sample(data[,3],n,replace=F)# d4<-sample(data[,4],n,replace=F)# shuffled_score<-sum((d1+d2+d3+d4)^2) # randomisations<-c(randomisations,shuffled_score)# } # hist(randomisations)# abline(v=raw_score,lwd=2,col=2)# text(528,12000,"observed",srt=270,col="red")# equals<-length(which(randomisations==raw_score))# below<-length(which(randomisationsraw_score))+0.5*equals# text(480,16000,paste("less than = ",below,sep=""))# text(560,16000,paste("greater than = ",above,sep=""),xpd=NA)# text(560,12000,paste("P = ",round(above/(above+below),4)),col="blue") day_in_year<-c(0.9023162,1.9475052,1.9111508,17.947441,17.61407,19.799334,33.986637, 33.615097,39.078434,49.624847,50.02184,58.024895,65.629875,65.61569,76.569275,81.67198,82.0228,95.510284,97.67119,98.09835,114.81702,113.7471,113.76201,133.79655,129.08539,145.89676,145.54774,152.80917,161.92287,171.76254,177.9617,191.43936,194.37025,209.6631,209.6791,228.60011,225.7012,247.55856,251.5441,254.8538,270.53964,270.5327,286.1364,286.54538,289.83145,302.56824,302.5624,308.39548,318.93423,318.5678,327.66803,334.57318,334.956,346.61084,350.22012,350.6073) # log_number<-c(1.434,0.454,-0.294,-0.123,0.513,0.506,0.109,-0.041,-0.056, -0.415, 0.258,0.034,-0.887,-1.179,-0.655,-0.596,-0.872,-0.678,-1.187, 0.108,-0.671,-0.199,0.107,0.099,0.601,1.73,2.044,1.55,1.69,1.782,1.916,1.909,2.253,2.118,2.447,2.012,2.327,2.35,1.907,2.55,3.01,2.86,1.63,2.55,2.707,2.445,2.325,2.422,1.906,1.86,1.726,1.397,1.779,1.741,1.0525,1.524) plot(day_in_year,log_number) # lines(lowess(day_in_year, log_number, f = .2, iter=5, delta=1),col="red",lwd=2) m<-lowess(day_in_year, log_number) # lines(m,col="green",lwd=2) model<-loess(log_number ~ day_in_year) # xv<-seq(0,366,0.1) # create a sequence 0:366 in increments of 0.1 # yv<-predict(model,xv) #use the loess model to predict y values for the xv date values # lines(xv,yv,col="blue",lwd=2) parasitoid<-c("Bracon1","Bracon1","Bracon1","Bracon1","Bracon1","agathidine1", "agathidine1","agathidine1","agathidine1","Cotesia1","Cotesia1","Cotesia1", "Cotesia1","Cotesia2","Cotesia3","Cotesia3","Cotesia3","Cotesia4","Aleiodes1", "Aleiodes2","Macrocentrus1") # host<-c("Nymphalidae2","Erebidae1","Erebidae1","Erebidae1","Erebidae1","Erebidae2", "Erebidae2","Noctuid1","Lycaenidae1","Hesperidae1","Hesperidae1","Hesperidae1", "Hesperidae2","Noctuidae2","Geometridae1","Geometridae2","Erebidae1", "Erebidae1","Erebidae1","Erebidae3", "Erebidae4") l<-levels(as.factor(parasitoid)) # m<-levels(as.factor(host)) plot(NULL, axes=F,xlab="",ylab="",xlim=c(0,length(parasitoid)),ylim=c(0,5)) P<-NULL # H<-NULL # for(i in 1:length(l)) P<-c(P,length(which(parasitoid==l[i]))) # for(i in 1:length(m)) H<-c(H,length(which(host==m[i]))) # S<-cumsum(P) # S1<-c(0,S) # T<-cumsum(H) # T1<-c(0,T) # for(i in 1:length(l)) rect(S1[i],3.9,S1[i+1],5.7,col="red",lwd=6,border="white") # for(j in 1:length(m)) rect(T1[j],0.7,T1[j+1],2,col="green",lwd=6,border="white") # midT<-T1[1:length(T)]+(T-T1[1:length(T)])/2 # midS<-S1[1:length(S)]+(S-S1[1:length(S)])/2 for(i in 1:length(l)) text(midS[i],4.1,l[i],srt=90,cex=0.8,xpd=NA,adj=0,font=3) # for(i in 1:length(m)) text(midT[i],1.9,m[i],srt=270,cex=0.8,xpd=NA,adj=0,font=3) # for(k in 1:length(l)){ # stepping through each different parasitoid # pp<-which(parasitoid==l[k]) # np<-length(pp) # hostsofparasitoid<-host[pp] # links<-length(levels(as.factor(hostsofparasitoid))) # for(q in 1:links){ # noRearings<-length(which(hostsofparasitoid==levels(as.factor(hostsofparasitoid))[q])) # temph<-which(m==levels(as.factor(hostsofparasitoid))[q]) # lines(c(midS[k],midT[temph]),c(3.8,2.2),col="grey50",lwd=noRearings*2,lend=3) # } # } reps<-100# generations<-500# start<-c(rep("a1",200),rep("a2",200),rep("a3",200),rep("a4",200),rep("a5",200))# allFreq<-start/1000# # a<-sample(start,1000)# # allele1<-list(a[1:500])# allele2<-list(a[501:1000])# genotypes<-rbind(allele1,allele2) genotypes a allele1 allele2 allele1<-c(a[1:500])# allele2<-c(a[501:1000])# genotypes<-rbind(allele1,allele2)# genotypes i<-1 j<-1 temp<-unlist(sample(genotypes[1:2,],25)) temp temp<-unlist(sample(genotypes[1:2,],50)) temp sample(genotypes[1:2,],50) n<-sample(1:500,25)# temp<-c(genotypes[1,n],genotypes[2,n])*20 genotypes[1, n] genotypes[2, n] temp<-c(genotypes[1,n],genotypes[2,n]) temp<-rep(c(genotypes[1,n],genotypes[2,n]),20) temp allele1<-c(temp[1:500])# allele2<-c(temp[501:1000])# genotypes<-rbind(allele1,allele2) levels(temp) levels(as.factor(temp)) reps<-100# generations<-500# start<-c(rep("a1",200),rep("a2",200),rep("a3",200),rep("a4",200),rep("a5",200))# # a<-sample(start,1000)# # allele1<-c(a[1:500])# allele2<-c(a[501:1000])# genotypes<-rbind(allele1,allele2)# genotypes# number_remaining<-matrix(NA,nrow=reps,ncol=generations)# # for(i in 1:reps){# for(j in 1:generations){# n<-sample(1:500,25)# temp<-sample(rep(c(genotypes[1,n],genotypes[2,n]),20),1000)# number_remaining[i,j]<-length(levels(as.factor(temp)))# allele1<-c(temp[1:500])# allele2<-c(temp[501:1000])# genotypes<-rbind(allele1,allele2)# }}# # plot(NULL,ylim=c(0,5),xlim=c(1,generations),xlab="Generation",ylab="Mean number of alleles in population")# # for(i in 1:generations) lines(c(i,i),c(mean(number_remaining[,i])-(1.96*sqrt(var(number_remaining[,i]))),mean(number_remaining[,i])+(1.96*sqrt(var(number_remaining[,i])))),pch=19,col="grey60")# for(i in 1:generations) points(i,mean(number_remaining[,i]),pch=16,col="blue",main="Mean alleles remaining (95% confidence intervals)") genotypes number_remaining[1,] number_remaining[2,] number_remaining[3,] i<-1 j<-1 start<-c(rep("a1",200),rep("a2",200),rep("a3",200),rep("a4",200),rep("a5",200))# # a<-sample(start,1000)# # allele1<-c(a[1:500])# allele2<-c(a[501:1000])# genotypes<-rbind(allele1,allele2)# genotypes number_remaining<-matrix(NA,nrow=reps,ncol=generations) n<-sample(1:500,25) n<-sample(1:500,25) n c(genotypes[1,n],genotypes[2,n]) c(genotypes[1,n],genotypes[1,n]) start<-c(rep("a1",200),rep("a2",200),rep("a3",200),rep("a4",200),rep("a5",200))# # a<-sample(start,1000)# # allele1<-c(a[1:500])# allele2<-c(a[501:1000])# genotypes<-rbind(allele1,allele2)# genotypes temp<-sample(rep(c(genotypes[1,n],genotypes[2,n]),20),1000) temp length(levels(as.factor(temp))) allele1<-c(temp[1:500])# allele2<-c(temp[501:1000])# genotypes<-rbind(allele1,allele2) genotypes number_remaining number_remaining[i,j]<-length(levels(as.factor(temp))) head(number_remaining)# ) head(number_remaining) j<-2 n<-sample(1:500,25)# temp<-sample(rep(c(genotypes[1,n],genotypes[2,n]),20),1000)# allele1<-c(temp[1:500])# allele2<-c(temp[501:1000])# genotypes<-rbind(allele1,allele2)# number_remaining[i,j]<-length(levels(as.factor(temp))) head(number_remaining,2) j<-3 n<-sample(1:500,25)# temp<-sample(rep(c(genotypes[1,n],genotypes[2,n]),20),1000)# allele1<-c(temp[1:500])# allele2<-c(temp[501:1000])# genotypes<-rbind(allele1,allele2)# number_remaining[i,j]<-length(levels(as.factor(temp))) head(number_remaining,2) start<-c(rep("a1",200),rep("a2",200),rep("a3",200),rep("a4",200),rep("a5",200))# # a<-sample(start,1000)# # allele1<-c(a[1:500])# allele2<-c(a[501:1000])# genotypes<-rbind(allele1,allele2)# genotypes# number_remaining<-matrix(NA,nrow=reps,ncol=generations)# # for(i in 1:reps){# for(j in 1:generations){# n<-sample(1:500,25)# temp<-sample(rep(c(genotypes[1,n],genotypes[2,n]),20),1000)# allele1<-c(temp[1:500])# allele2<-c(temp[501:1000])# genotypes<-rbind(allele1,allele2)# number_remaining[i,j]<-length(levels(as.factor(temp)))# }}# # plot(NULL,ylim=c(0,5),xlim=c(1,generations),xlab="Generation",ylab="Mean number of alleles in population")# # for(i in 1:generations) lines(c(i,i),c(mean(number_remaining[,i])-(1.96*sqrt(var(number_remaining[,i]))),mean(number_remaining[,i])+(1.96*sqrt(var(number_remaining[,i])))),pch=19,col="grey60")# for(i in 1:generations) points(i,mean(number_remaining[,i]),pch=16,col="blue",main="Mean alleles remaining (95% confidence intervals)") number_remaining<-matrix(NA,nrow=reps,ncol=generations)# # for(i in 1:reps){# start<-c(rep("a1",200),rep("a2",200),rep("a3",200),rep("a4",200),rep("a5",200))# a<-sample(start,1000)# allele1<-c(a[1:500])# allele2<-c(a[501:1000])# genotypes<-rbind(allele1,allele2)# # for(j in 1:generations){# n<-sample(1:500,25)# temp<-sample(rep(c(genotypes[1,n],genotypes[2,n]),20),1000)# allele1<-c(temp[1:500])# allele2<-c(temp[501:1000])# genotypes<-rbind(allele1,allele2)# number_remaining[i,j]<-length(levels(as.factor(temp)))# }}# # plot(NULL,ylim=c(0,5),xlim=c(1,generations),xlab="Generation",ylab="Mean number of alleles in population")# # for(i in 1:generations) lines(c(i,i),c(mean(number_remaining[,i])-(1.96*sqrt(var(number_remaining[,i]))),mean(number_remaining[,i])+(1.96*sqrt(var(number_remaining[,i])))),pch=19,col="grey60")# for(i in 1:generations) points(i,mean(number_remaining[,i]),pch=16,col="blue",main="Mean alleles remaining (95% confidence intervals)") number_remaining<-matrix(NA,nrow=reps,ncol=generations)# # for(i in 1:reps){# start<-c(rep("a1",200),rep("a2",200),rep("a3",200),rep("a4",200),rep("a5",200))# a<-sample(start,1000)# allele1<-c(a[1:500])# allele2<-c(a[501:1000])# genotypes<-rbind(allele1,allele2)# # for(j in 1:generations){# n<-sample(1:500,25)# temp<-sample(rep(c(genotypes[1,n],genotypes[2,n]),20),1000)# allele1<-c(temp[1:500])# allele2<-c(temp[501:1000])# genotypes<-rbind(allele1,allele2)# number_remaining[i,j]<-length(levels(as.factor(temp)))# }}# # plot(NULL,ylim=c(0,5),xlim=c(1,generations),xlab="Generation",ylab="Mean number of alleles in population")# # for(i in 1:generations) lines(c(i,i),c(mean(number_remaining[,i])-(1.96*sqrt(var(number_remaining[,i]))),mean(number_remaining[,i])+(1.96*sqrt(var(number_remaining[,i])))),pch=19,col="grey60")# for(i in 1:generations) points(i,mean(number_remaining[,i]),pch=16,col="blue",main="Mean alleles remaining (95% confidence intervals)")