one.year.herblow<-function(Ind,YQ1,YQ2,YQ3,Infest,PropD,PropDorm){ Y1<-YQ1 Y2<-YQ2 Y3<-YQ3 Ind<-numeric(101) Ind<-Ind N<-1 I<-Infest PD<-PropD PDorm<-PropDorm P.d<-numeric(N) Dbankmort<-numeric(N) Ibankmort<-numeric(N) Rbankmort<-numeric(N) OYieldd<-numeric(N) OYieldi<-numeric(N) OYieldr<-numeric(N) discsurv<-numeric(N) intersurv<-numeric(N) raysurv<-numeric(N) dormbankbute<-numeric(N) raybankbute<-numeric(N) discemergence<-numeric(N) rayemergeyr1<-numeric(N) interemergeyr1<-numeric(N) interemergeyr2<-numeric(N) numdisc<-capittotal1*PD numdormant<-capittotal1*PDorm Yielda<-cbind(Disc=round(numdisc,0),Dorm=round(numdormant,0)) numinf<-sample(1:30,Infest,replace=FALSE) ############ ASSIGN DAMAGE PER INFESTATION for (i in 1:length(numinf)){ Yielda[numinf[i],1] <- Yielda[numinf[i],1]-(Yielda[numinf[i],1] * rnorm(1,0.535,0)) #(Non-dormant D4-D14)# Yielda[numinf[i],2] <- Yielda[numinf[i],2]-(Yielda[numinf[i],2] * rnorm(1,0.18,0)) #(Dormant D4-D14)# if (Yielda[numinf[i],1]<0){ Yielda[numinf[i],1]<-0 } if (Yielda[numinf[i],2]<0){ Yielda[numinf[i],2]<-0 } } Yieldd1<-Yielda[,1] Yielddorm1<-Yielda[,2] ###### MORTALITY ######## Discbankmort<-0.294 Dormbankmort<-0.764 Discbankmort<-abs(Discbankmort) if(Ibankmort>1){Ibankmort<-1} ###### GERMINATION OUT OF VIABLE ######## Germ.d<-1 devGerm.d<-0 Germ.dorm1<-0.1183 Germ.dorm2<-0.932 Germ.dorm3<-1 OYieldd<-sum(Yieldd1) OYieldi<-sum(Yielddorm1) DiscViable<-OYieldd*(1-Discbankmort) DormViable1<-OYieldi*(1-Dormbankmort) dormemergeyr1<-sum(Germ.dorm1*DormViable1) discemergence<-sum(DiscViable) DormViable2<-DormViable1-dormemergeyr1 dormemergeyr2<-sum(DormViable2*Germ.dorm2) dormemergeyr3<-DormViable1-(dormemergeyr2+dormemergeyr1) Sd<-(discemergence) Si<-(dormemergeyr1) Ri<-(dormemergeyr2) Ti<-sum(dormemergeyr3) ### end changes ### S<-sum(Sd,Si) R<-sum(Ri) ### what is P? year one survival! either 0 or 1 0= death of all that germinated that year ### ### what is Q? year two (next year) either 0 or 1 1=survival of all that germinated that year ### S<-S*Y1 R<-R*Y2 T<-Ti*Y3 Lambda<-(S[1]+R[1]+T[1]) return(Lambda) } one.year.herbhigh<-function(Ind,YQ1,YQ2,YQ3,Infest,PropD,PropDorm){ Y1<-YQ1 Y2<-YQ2 Y3<-YQ3 Ind<-Ind N<-1 I<-Infest PD<-PropD PDorm<-PropDorm P.d<-numeric(N) Dbankmort<-numeric(N) Ibankmort<-numeric(N) Rbankmort<-numeric(N) OYieldd<-numeric(N) OYieldi<-numeric(N) OYieldr<-numeric(N) discsurv<-numeric(N) intersurv<-numeric(N) raysurv<-numeric(N) dormbankbute<-numeric(N) raybankbute<-numeric(N) discemergence<-numeric(N) rayemergeyr1<-numeric(N) interemergeyr1<-numeric(N) interemergeyr2<-numeric(N) numdisc<-capittotal1*PD numdormant<-capittotal1*PDorm Yielda<-cbind(Disc=round(numdisc,0),Dorm=round(numdormant,0)) numinf<-sample(1:30,Infest,replace=FALSE) ############ ASSIGN DAMAGE PER INFESTATION for (i in 1:length(numinf)){ Yielda[numinf[i],1] <- Yielda[numinf[i],1]-(Yielda[numinf[i],1] * rnorm(1,0.7980,0)) #(Non-dormant D4 only)# Yielda[numinf[i],2] <- Yielda[numinf[i],2]-(Yielda[numinf[i],2] * rnorm(1,0.358,0)) #(Dormant D4 only)# if (Yielda[numinf[i],1]<0){ Yielotda[numinf[i],1]<-0 } if (Yielda[numinf[i],2]<0){ Yielda[numinf[i],2]<-0 } } Yieldd1<-Yielda[,1] Yielddorm1<-Yielda[,2] ######MORTALITY######## Discbankmort<-0.294 Dormbankmort<-0.764 Discbankmort<-abs(Discbankmort) if(Ibankmort>1){Ibankmort<-1} ###### GERMINATION OUT OF VIABLE ######## Germ.d<-1 devGerm.d<-0 Germ.dorm1<-0.1183 Germ.dorm2<-0.932 Germ.dorm3<-1 OYieldd<-sum(Yieldd1) OYieldi<-sum(Yielddorm1) DiscViable<-OYieldd*(1-Discbankmort) DormViable1<-OYieldi*(1-Dormbankmort) dormemergeyr1<-sum(Germ.dorm1*DormViable1) discemergence<-sum(DiscViable) DormViable2<-DormViable1-dormemergeyr1 dormemergeyr2<-sum(DormViable2*Germ.dorm2) dormemergeyr3<-DormViable1-(dormemergeyr2+dormemergeyr1) Sd<-(discemergence) Si<-(dormemergeyr1) Ri<-(dormemergeyr2) Ti<-sum(dormemergeyr3) ### end changes ### S<-sum(Sd,Si) R<-sum(Ri) ### what is P? year one survival! either 0 or 1 0= death of all that germinated that year ### ### what is Q? year two (next year) either 0 or 1 1=survival of all that germinated that year ### S<-S*Y1 R<-R*Y2 T<-Ti*Y3 Lambda<-(S[1]+R[1]+T[1]) return(Lambda) } ######## begin with 0 Disc and only intermediate and rays. ###### ######## end with only Disc #### Pred<-read.table("Achene.cor.predict",header=TRUE) Pred$Total<-Pred$Disc+Pred$Inter+Pred$Ray predtotal<-lm(Total~Florets,data=Pred) ncapit1<-round(rnorm(1,30,0),0) ncapit21<-numeric(sum(ncapit1)) ncap1<-length(ncapit21) capitfloret1<-numeric(ncap1) capitfloret1<-round(rnorm(ncap1,40,0),0) newfloret1<-data.frame(Florets=capitfloret1) capittotal<-predict(predtotal,interval="predict",newdata=newfloret1) #capittotal1<-data.frame(Total=(capittotal[,1]+runif(length(capittotal[,1]),-87,87)))# capittotal1<-data.frame(Total=(capittotal[,1])) REVISEDEnviro.geomean<-matrix(nrow=101,ncol=2) REVISEDEnviro.sd<-matrix(nrow=101,ncol=2) REVISEDEnviro.arithmean<-matrix(nrow=101,ncol=2) SDREVISEDEnviro.geomean<-matrix(nrow=101,ncol=2) SDREVISEDEnviro.sd<-matrix(nrow=101,ncol=2) SDREVISEDEnviro.arithmean<-matrix(nrow=101,ncol=2) P<-seq(0,1,.01) Q<-1-P Infest<-30 YQ0<-rep(1,102) for (e in 1:101){ REVISEDFitness<-matrix(nrow=100,ncol=4) REVISEDLambda<-matrix(nrow=100,ncol=101) REVISEDSd<-matrix(nrow=100,ncol=4) ###### replicates go here ##### for (j in 3:102){ for (i in 1:101){ ifelse(BG[j-2,e]==0, REVISEDLambda[j-2,i]<-one.year.herbhigh(1,YQ0[j-2],YQ0[j-1],YQ0[j],Infest,P[i],Q[i]), REVISEDLambda[j-2,i]<-one.year.herblow(1,YQ0[j-2],YQ0[j-1],YQ0[j],Infest,P[i],Q[i])) } } #REVISEDLambda<-REVISEDLambda[-1,] # REVISEDLambda.ARITHmean<-numeric(101) REVISEDLambda.GEOmean<-numeric(101) REVISEDLambda.sd<-numeric(101) for (i in 1:101){ REVISEDLambda.ARITHmean[i]<-mean(REVISEDLambda[,i]) REVISEDLambda.sd[i]<-sd(REVISEDLambda[,i]) REVISEDLambda.GEOmean[i]<-exp(mean(log(REVISEDLambda[-1,i]))) } o<-cbind("Geo"=REVISEDLambda.GEOmean,"Arith"=REVISEDLambda.ARITHmean,"SD"=REVISEDLambda.sd,"Disc"=P) REVISEDFitness<-o[o[,1] == max(o[,1]), ] REVISEDSd<-o[o[,3] == min(o[,3]), ] REVISEDEnviro.geomean[e,]<-cbind(mean(REVISEDFitness[1]),mean(REVISEDFitness[4])) REVISEDEnviro.arithmean[e,]<-cbind(mean(REVISEDFitness[2]),mean(REVISEDFitness[4])) REVISEDEnviro.sd[e,]<-cbind(mean(REVISEDFitness[3]),mean(REVISEDFitness[4])) SDREVISEDEnviro.geomean[e,]<-cbind(mean(REVISEDSd[1]),mean(REVISEDSd[4])) SDREVISEDEnviro.arithmean[e,]<-cbind(mean(REVISEDSd[2]),mean(REVISEDSd[4])) SDREVISEDEnviro.sd[e,]<-cbind(mean(REVISEDSd[3]),mean(REVISEDSd[4])) } NEWEnviro.arithmean.DormHerb0.00<-REVISEDEnviro.arithmean NEWEnviro.geomean.DormHerb0.00<-REVISEDEnviro.geomean NEWEnviro.sd.DormHerb0.00<-REVISEDEnviro.sd SDNEWEnviro.arithmean.DormHerb0.00<-SDREVISEDEnviro.arithmean SDNEWEnviro.geomean.DormHerb0.00<-SDREVISEDEnviro.geomean SDNEWEnviro.sd.DormHerb0.00<-SDREVISEDEnviro.sd save("NEWEnviro.arithmean.DormHerb0.00","NEWEnviro.geomean.DormHerb0.00","NEWEnviro.sd.DormHerb0.00",file="NEW.DormHerb0.00.3year") save("SDNEWEnviro.arithmean.DormHerb0.00","SDNEWEnviro.geomean.DormHerb0.00","SDNEWEnviro.sd.DormHerb0.00",file="SDNEW.DormHerb0.00.3year")