############################################### # Normal and skew normal across parameter values (Tsds and CTsds) ############################################### erf <- function(x) 2 * pnorm(x * sqrt(2)) - 1 require(Hmisc) # for wtd.var require(pracma) # for erf ntemp <- 400 # number of temperatures (including Tref, which we add manually) Tref<-15 Ts <- seq(0,40, length.out=ntemp) #ntemp<-3 #Ts<-c(14.9,15,15.1) if (!Tref %in% Ts){ Ts <- append(Ts,Tref) #We make sure Tref is in the list ntemp<-ntemp+1 print("We need to add Tref manually") } Ts<-sort(Ts,decreasing=FALSE) nspp <- 1000 nreps <- 1000 # how many repeats at each combination of Tsds and CTsds Tsds <- seq(1,15,by=2) # set of species thermal standard deviations##for this mixed case, it now sets the width of the niche for species in the pool CTsds <- seq(0,30, by=0.5) # set of species pool thermal standard deviations lambda <- -10 STI_S <- array(NA, dim=c(length(nspp))) CTIN_T <- array(NA, dim=c(length(Ts), length(Tsds), length(CTsds))) CTIS_T <- array(NA, dim=c(length(Ts), length(Tsds), length(CTsds))) CTDivsN_T <- array(NA, dim=c(length(Ts), length(Tsds), length(CTsds))) CTDivsS_T <- array(NA, dim=c(length(Ts), length(Tsds), length(CTsds))) sCTIsN_T <- array(NA, dim=c(length(Ts), length(Tsds), length(CTsds))) sCTIsS_T <- array(NA, dim=c(length(Ts), length(Tsds), length(CTsds))) ShN_T <- array(NA, dim=c(length(Ts), length(Tsds), length(CTsds))) ShS_T <- array(NA, dim=c(length(Ts), length(Tsds), length(CTsds))) sCTIsN <- array(NA, dim=c(length(Tsds), length(CTsds))) sCTIsS <- array(NA, dim=c(length(Tsds), length(CTsds))) CTDivsN <- array(NA, dim=c(length(Tsds), length(CTsds))) CTDivsS <- array(NA, dim=c(length(Tsds), length(CTsds))) sCTIpredN <- array(0, dim=c(length(Tsds), length(CTsds))) sCTIpredS <- array(0, dim=c(length(Tsds), length(CTsds))) CTINmean_T<-array(0, dim=c(length(Ts), length(Tsds), length(CTsds))) CTISmean_T<-array(0, dim=c(length(Ts), length(Tsds), length(CTsds))) CTDivsNmean_T<-array(0, dim=c(length(Ts), length(Tsds), length(CTsds))) CTDivsSmean_T<-array(0, dim=c(length(Ts), length(Tsds), length(CTsds))) sCTIsNmean_T<-array(0, dim=c(length(Ts), length(Tsds), length(CTsds))) sCTIsSmean_T<-array(0, dim=c(length(Ts), length(Tsds), length(CTsds))) ShNmean_T<-array(0, dim=c(length(Ts), length(Tsds), length(CTsds))) ShSmean_T<-array(0, dim=c(length(Ts), length(Tsds), length(CTsds))) set.seed(642)#Not strictly needed to generate the random numbers, but ensures reproducibility. for(i in 1:length(Tsds)){ # for all species thermal standard deviations print(paste(i, 'of', length(Tsds))) for(j in 1:length(CTsds)){ #and for all species pool thermal standard deviations print(paste(' ', j, 'of', length(CTsds))) for(r in 1:nreps){ # for all replicates of each combination of standard deviation and pool standard deviation Topts <- rnorm(nspp, mean=Tref, sd=CTsds[j]) # Toptimum for each species, same a species temperature index for normal # STI_S <- Topts + Tsds[i]*(lambda/sqrt(1+lambda^2))*sqrt(2/pi) # species temperature index for skew normal (the mean) tresponseN <- matrix(NA, ncol=ntemp, nrow=nspp) # normal tresponseS <- matrix(NA, ncol=ntemp, nrow=nspp) # skew normal sigma_vector<-runif(nspp,0,Tsds[i]) for(k in 1:nspp){ sigma<-sigma_vector[k] STI_S[k] <- Topts[k] + sigma*(lambda/sqrt(1+lambda^2))*sqrt(2/pi) tresponseN[k,] <- (1/sqrt(2*pi))*(1/sigma)*exp(-(Ts-Topts[k])^2/(2*sigma^2))*(1+erf(0.*(Ts-Topts[k])/(sqrt(2)*sigma))) tresponseS[k,] <- (1/sqrt(2*pi))*(1/sigma)*exp(-(Ts-Topts[k])^2/(2*sigma^2))*(1+erf(lambda*(Ts-Topts[k])/(sqrt(2)*sigma))) } CTIN <- as.numeric(Topts %*% tresponseN)/colSums(tresponseN) # calculate CTI for normal CTIS <- as.numeric(STI_S %*% tresponseS)/colSums(tresponseS) # calculate CTI for skew normal using mean of skew normal for(k in 2:length(Ts)){ inds <- c(k-1,k) sCTIsN_T[k,i,j]<-diff(as.numeric(CTIN[inds]))/diff(Ts[inds]) # calculate CTI sensitivity for normal case sCTIsS_T[k,i,j]<-diff(as.numeric(CTIS[inds]))/diff(Ts[inds]) # calculate CTI sensitivity for skew normal case } for(k in 1:length(Ts)){ CTIN_T[k,i,j]<-CTIN[k] CTIS_T[k,i,j]<-CTIS[k] CTDivsN_T[k,i,j]<-sqrt(((Topts-CTIN[k])^2 %*% tresponseN[,k])/sum(tresponseN[,k])) # for normal CTDivsS_T[k,i,j]<-sqrt(((STI_S-CTIS[k])^2 %*% tresponseS[,k])/sum(tresponseS[,k])) # using mean of skew normal temp<-tresponseN[which(tresponseN[,k]>0),k]/sum(tresponseN[which(tresponseN[,k]>0),k]) #in case there's any zero temp<-temp[which(temp>0)] #in case the element is so small that is a numerical zero ShN_T[k,i,j] <- sum(-temp*log10(temp)) # calculate Shannon biodiversity index for normal case temp<-tresponseS[which(tresponseS[,k]>0),k]/sum(tresponseS[which(tresponseS[,k]>0),k]) #in case there's any zero temp<-temp[which(temp>0)] #in case the element is so small that is a numerical zero ShS_T[k,i,j] <- sum(-temp*log10(temp)) # calculate Shannon biodiversity index for skew normal case # adding the information of all observables across replicates for the calculation of the mean later: CTINmean_T[k,i,j] <- CTINmean_T[k,i,j] + CTIN_T[k,i,j] CTISmean_T[k,i,j] <- CTISmean_T[k,i,j] + CTIS_T[k,i,j] CTDivsNmean_T[k,i,j] <- CTDivsNmean_T[k,i,j] + CTDivsN_T[k,i,j] CTDivsSmean_T[k,i,j] <- CTDivsSmean_T[k,i,j] + CTDivsS_T[k,i,j] sCTIsNmean_T[k,i,j] <- sCTIsNmean_T[k,i,j] + sCTIsN_T[k,i,j] sCTIsSmean_T[k,i,j] <- sCTIsSmean_T[k,i,j] + sCTIsS_T[k,i,j] ShNmean_T[k,i,j] <- ShNmean_T[k,i,j] + ShN_T[k,i,j] ShSmean_T[k,i,j] <- ShSmean_T[k,i,j] + ShS_T[k,i,j] } n<-which(Ts==Tref) # specific case in which climate matches the chosen Tref value smean<-median(sigma_vector) # we choose the median from all possible sigmas sCTIpredN[i,j] <- sCTIpredN[i,j] + (CTDivsN_T[n,i,j]^2/smean^2)+(0./smean)*sqrt(2/pi)*(Tref+sqrt(2/pi)*smean*0./sqrt(1+0.^2)-CTIN_T[n,i,j])*(sqrt(smean^2+CTsds[j]^2)/sqrt(smean^2+(1+0.^2)*CTsds[j]^2)) sCTIpredS[i,j] <- sCTIpredS[i,j] + (CTDivsS_T[n,i,j]^2/smean^2)+(lambda/smean)*sqrt(2/pi)*(Tref+sqrt(2/pi)*smean*lambda/sqrt(1+lambda^2)-CTIS_T[n,i,j])*(sqrt(smean^2+CTsds[j]^2)/sqrt(smean^2+(1+lambda^2)*CTsds[j]^2)) } } } # summarize across replicates CTINmean_T <- CTINmean_T/nreps CTISmean_T <- CTISmean_T/nreps CTDivsNmean_T <- CTDivsNmean_T/nreps CTDivsSmean_T <- CTDivsSmean_T/nreps sCTIsNmean_T <- sCTIsNmean_T/nreps sCTIsSmean_T <- sCTIsSmean_T/nreps ShNmean_T <- ShNmean_T/nreps ShSmean_T <- ShSmean_T/nreps n<-which(Ts==Tref) sCTIsNmean <- sCTIsNmean_T[n,,] sCTIsSmean <- sCTIsSmean_T[n,,] sCTIpredNmean <- sCTIpredN/nreps sCTIpredSmean <- sCTIpredS/nreps CTDivsNmean <- CTDivsNmean_T[n,,] CTDivsSmean <- CTDivsSmean_T[n,,] ShNmean <- ShNmean_T[n,,] ShSmean <- ShSmean_T[n,,] ##Export to files: for(i in 1:length(Tsds)){ # which Tsds to choose sigma=sprintf("%05.2f",Tsds[i]) print(paste0(sigma," ",i)) name=paste0("slope_vs_CTIDiv_sigma_",sigma,".txt") OUT=cbind(CTDivsNmean[i,], sCTIsNmean[i,],sCTIpredNmean[i,],ShNmean[i,],CTDivsSmean[i,], sCTIsSmean[i,],sCTIpredSmean[i,],ShSmean[i,],Tref,Tsds[i],CTsds) write.table(OUT, name, sep="\t", row.names=FALSE, col.names=FALSE) for(j in 1:length(CTsds)){ breadthT=sprintf("%05.2f",CTsds[j]) print(paste0(sigma," ",i," ",breadthT," ",j)) name=paste0("all_curves_vs_T_sigma_",sigma,"_breadthT_",breadthT,".txt") OUT=cbind(Ts,CTINmean_T[,i,j],CTDivsNmean_T[,i,j],sCTIsNmean_T[,i,j],ShNmean_T[,i,j],CTISmean_T[,i,j],CTDivsSmean_T[,i,j],sCTIsSmean_T[,i,j],ShSmean_T[,i,j],Tref,Tsds[i],CTsds[j]) write.table(OUT, name, sep="\t", row.names=FALSE, col.names=FALSE) } }