Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Francomano_et_al_2020_Temporal_Variability/Question_2_Analysis.R
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
604 lines (566 sloc)
37.9 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ###################################################################### | |
| # Code from Francomano et al. (Ecological Indicators 2020) # | |
| ###################################################################### | |
| # Question 2 Analysis # | |
| ###################################################################### | |
| # Created by: Dante Francomano (dfrancom@purdue.edu) | |
| # Date: 20180518 | |
| ###################################################################### | |
| ###################################################################### | |
| # Set working directory and load packages | |
| setwd("/Users/dantefrancomano/Documents/Purdue/Projects/Temporal_Variability_Duty_Cycle/Datasets") | |
| library(dplyr) | |
| library(lme4) | |
| library(MuMIn) | |
| library(lattice) | |
| #################################################################################################### | |
| # Create dataframe for results | |
| #################################################################################################### | |
| sites<-read.csv("/Users/dantefrancomano/Documents/Purdue/Projects/Temporal_Variability_Duty_Cycle/Metadata.csv")[,2] | |
| timeOfDay<-c("dawn","daytime","dusk","nighttime") | |
| index<-c("ACI","BI","Hf") | |
| subdivisions<-function(x,y) | |
| { | |
| possibleFactors<-seq(1,x) | |
| nFactors<-length(which((x/possibleFactors)==trunc(x/possibleFactors)&(y/possibleFactors)==trunc(y/possibleFactors)&(y/x)==trunc(y/x))) | |
| nFactors[which(nFactors==0)]<-1 | |
| if(x==1|x>(y/2)) return(1) | |
| if(x<=(y/2)) return(nFactors) | |
| } | |
| factors<-function(x) | |
| { | |
| possibleFactors<-seq(1,x) | |
| which((x/possibleFactors)==trunc(x/possibleFactors)) | |
| } | |
| numberOfSubdivisions12<-sapply(1:12,subdivisions,12) | |
| numberOfSubdivisions60<-sapply(1:60,subdivisions,60) | |
| totalAnalysisTime12<-as.numeric(unlist(apply(data.frame(1:12,numberOfSubdivisions12),1,function(y) rep(y[1],times=1,each=y[2])))) | |
| totalAnalysisTime60<-as.numeric(unlist(apply(data.frame(1:60,numberOfSubdivisions60),1,function(y) rep(y[1],times=1,each=y[2])))) | |
| specificSubdivisions12<-unlist(apply(data.frame(1:12,numberOfSubdivisions12),1,function(x) ifelse(x[2]==1,return(1),return(factors(x[1]))))) | |
| specificSubdivisions60<-unlist(apply(data.frame(1:60,numberOfSubdivisions60),1,function(x) ifelse(x[2]==1,return(1),return(factors(x[1]))))) | |
| results12<-data.frame(Location=rep(sites,each=length(totalAnalysisTime12)*4), | |
| Time_of_Day=rep(timeOfDay,times=8,each=length(totalAnalysisTime12)), | |
| Total_Analysis_Time=rep(totalAnalysisTime12,times=32), | |
| Number_of_Subdivisions=rep(specificSubdivisions12,times=32), | |
| ACI_Deviation=NA,BI_Deviation=NA,Hf_Deviation=NA) | |
| results60<-data.frame(Location=rep(sites,each=length(totalAnalysisTime60)*4), | |
| Time_of_Day=rep(timeOfDay,times=8,each=length(totalAnalysisTime60)), | |
| Total_Analysis_Time=rep(totalAnalysisTime60,times=32), | |
| Number_of_Subdivisions=rep(specificSubdivisions60,times=32), | |
| ACI_Deviation=NA,BI_Deviation=NA,Hf_Deviation=NA) | |
| # Read in data | |
| datasetList<-list.files() | |
| datasets<-lapply(datasetList,function(x) | |
| { | |
| read.csv(x) | |
| }) | |
| names(datasets)<-substr(datasetList,1,(nchar(datasetList)-4)) | |
| #################################################################################################### | |
| # Calculate global index means for data used in analysis | |
| #################################################################################################### | |
| oneMinuteTimeOfDayValues<-NULL | |
| dawnDatasets<-datasets[which(grepl(paste("Dawn","_1-minute_",sep=""),names(datasets)))] | |
| daytimeDatasets<-datasets[which(grepl(paste("Daytime","_1-minute_2-hour",sep=""),names(datasets)))] | |
| duskDatasets<-datasets[which(grepl(paste("Dusk","_1-minute_",sep=""),names(datasets)))] | |
| nighttimeDatasets<-NULL | |
| for(i in 1:8) nighttimeDatasets[i]<-list(datasets[which(grepl(paste("Nighttime","_1-minute_2-hour_",sites[i],sep=""),names(datasets)))][[1]][,-4]) | |
| timeOfDayDatasets<-c(dawnDatasets,daytimeDatasets,duskDatasets,nighttimeDatasets) | |
| for(i in 1:32) oneMinuteTimeOfDayValues<-rbind(oneMinuteTimeOfDayValues,timeOfDayDatasets[i][[1]]) | |
| colMeans(oneMinuteTimeOfDayValues[,4:6],na.rm=TRUE) | |
| #################################################################################################### | |
| # Calculate subsample error values | |
| #################################################################################################### | |
| # Begin site loop | |
| for(g in 1:length(sites)) | |
| { | |
| # Begin index loop | |
| for(h in 1:length(index)) | |
| { | |
| # Begin time of day loop | |
| for(i in 1:length(timeOfDay)) | |
| { | |
| # Create object for time of day- and value length-specific datasets | |
| if(i<4) | |
| { | |
| oneMinute<-datasets[which(grepl(paste("D",substr(timeOfDay[i],2,nchar(timeOfDay[i])),"_1-minute_",sites[g],sep=""),names(datasets)))][[1]] | |
| }else{ | |
| oneMinute<-datasets[which(grepl(paste("Nighttime_1-minute_2-hour_",sites[g],sep=""),names(datasets)))][[1]] | |
| oneMinute<-oneMinute[,-4] | |
| } | |
| # Conduct consecutive sampling measurements | |
| series12<-NULL | |
| for(j in 1:(dim(oneMinute)[1]/12)) | |
| { | |
| if(!is.na(sum(oneMinute[(j*12-11):(j*12),(3+h)]))) | |
| { | |
| mean12<-sum(oneMinute[(j*12-11):(j*12),(3+h)],na.rm=TRUE)/12 | |
| }else{ | |
| mean12<-sum(oneMinute[(j*12-11):(j*12),(3+h)],na.rm=TRUE)/11 | |
| } | |
| for(k in 1:12) | |
| { | |
| if(!is.na(sum(oneMinute[(j*12-11):(j*12-11+k-1),(3+h)]))) | |
| { | |
| cumMean<-sum(oneMinute[(j*12-11):(j*12-11+k-1),(3+h)],na.rm=TRUE)/k | |
| }else{ | |
| cumMean<-sum(oneMinute[(j*12-11):(j*12-11+k-1),(3+h)],na.rm=TRUE)/(k-1) | |
| } | |
| series12[(j-1)*12+k]<-abs(cumMean-mean12) | |
| } | |
| } | |
| consecMeans12<-(colMeans(matrix(series12,ncol=12,byrow=TRUE),na.rm=TRUE)) | |
| results12[which(results12[,1]==sites[g]&results12[,2]==timeOfDay[i]&results12[,4]==1),(4+h)]<-consecMeans12 | |
| series60<-NULL | |
| for(j in 1:(dim(oneMinute)[1]/60)) | |
| { | |
| if(!is.na(sum(oneMinute[(j*60-59):(j*60),(3+h)]))) | |
| { | |
| mean60<-sum(oneMinute[(j*60-59):(j*60),(3+h)],na.rm=TRUE)/60 | |
| }else{ | |
| mean60<-sum(oneMinute[(j*60-59):(j*60),(3+h)],na.rm=TRUE)/59 | |
| } | |
| for(k in 1:60) | |
| { | |
| if(!is.na(sum(oneMinute[(j*60-59):(j*60-59+k-1),(3+h)]))) | |
| { | |
| cumMean<-sum(oneMinute[(j*60-59):(j*60-59+k-1),(3+h)],na.rm=TRUE)/k | |
| }else{ | |
| cumMean<-sum(oneMinute[(j*60-59):(j*60-59+k-1),(3+h)],na.rm=TRUE)/(k-1) | |
| } | |
| series60[(j-1)*60+k]<-abs(cumMean-mean60) | |
| } | |
| } | |
| consecMeans60<-(colMeans(matrix(series60,ncol=60,byrow=TRUE),na.rm=TRUE)) | |
| results60[which(results60[,1]==sites[g]&results60[,2]==timeOfDay[i]&results60[,4]==1),(4+h)]<-consecMeans60 | |
| # Conduct distributed sampling measurements | |
| distributedResults<-results12[which(results12[,1]==sites[g]&results12[,2]==timeOfDay[i]&results12[,4]>1),] | |
| for(j in 1:dim(distributedResults)[1]) | |
| { | |
| tempResults<-NULL | |
| totalAnalysisTime<-distributedResults[j,3] | |
| numberOfSubdivisions<-distributedResults[j,4] | |
| sampleLength<-totalAnalysisTime/numberOfSubdivisions | |
| sampleInterval<-12/numberOfSubdivisions | |
| sampleContents<-NULL | |
| for(k in 1:numberOfSubdivisions) | |
| { | |
| consecutiveSequence<-seq(((k-1)*sampleInterval+1),((k-1)*sampleInterval+sampleLength)) | |
| sampleContents<-c(sampleContents,consecutiveSequence) | |
| } | |
| for(k in 1:(dim(oneMinute)[1]/12)) | |
| { | |
| if(!is.na(sum(oneMinute[(k*12-11):(k*12),(3+h)]))) | |
| { | |
| mean12<-sum(oneMinute[(k*12-11):(k*12),(3+h)],na.rm=TRUE)/12 | |
| }else{ | |
| mean12<-sum(oneMinute[(k*12-11):(k*12),(3+h)],na.rm=TRUE)/11 | |
| } | |
| if(!is.na(sum(oneMinute[((k-1)*12+sampleContents),(3+h)]))) | |
| { | |
| sampleMean<-sum(oneMinute[((k-1)*12+sampleContents),(3+h)],na.rm=TRUE)/totalAnalysisTime | |
| }else{ | |
| sampleMean<-sum(oneMinute[((k-1)*12+sampleContents),(3+h)],na.rm=TRUE)/(totalAnalysisTime-1) | |
| } | |
| tempResults[k]<-abs(sampleMean-mean12) | |
| } | |
| results12[which(results12[,1]==sites[g]&results12[,2]==timeOfDay[i]&results12[,3]==distributedResults[j,3]&results12[,4]==distributedResults[j,4]),(4+h)]<-mean(tempResults) | |
| } | |
| distributedResults<-results60[which(results60[,1]==sites[g]&results60[,2]==timeOfDay[i]&results60[,4]>1),] | |
| for(j in 1:dim(distributedResults)[1]) | |
| { | |
| tempResults<-NULL | |
| totalAnalysisTime<-distributedResults[j,3] | |
| numberOfSubdivisions<-distributedResults[j,4] | |
| sampleLength<-totalAnalysisTime/numberOfSubdivisions | |
| sampleInterval<-60/numberOfSubdivisions | |
| sampleContents<-NULL | |
| for(k in 1:numberOfSubdivisions) | |
| { | |
| consecutiveSequence<-seq(((k-1)*sampleInterval+1),((k-1)*sampleInterval+sampleLength)) | |
| sampleContents<-c(sampleContents,consecutiveSequence) | |
| } | |
| for(k in 1:(dim(oneMinute)[1]/60)) | |
| { | |
| if(!is.na(sum(oneMinute[(k*60-59):(k*60),(3+h)]))) | |
| { | |
| mean60<-sum(oneMinute[(k*60-59):(k*60),(3+h)],na.rm=TRUE)/60 | |
| }else{ | |
| mean60<-sum(oneMinute[(k*60-59):(k*60),(3+h)],na.rm=TRUE)/59 | |
| } | |
| if(!is.na(sum(oneMinute[((k-1)*60+sampleContents),(3+h)]))) | |
| { | |
| sampleMean<-sum(oneMinute[((k-1)*60+sampleContents),(3+h)],na.rm=TRUE)/totalAnalysisTime | |
| }else{ | |
| sampleMean<-sum(oneMinute[((k-1)*60+sampleContents),(3+h)],na.rm=TRUE)/(totalAnalysisTime-1) | |
| } | |
| tempResults[k]<-abs(sampleMean-mean60) | |
| } | |
| results60[which(results60[,1]==sites[g]&results60[,2]==timeOfDay[i]&results60[,3]==distributedResults[j,3]&results60[,4]==distributedResults[j,4]),(4+h)]<-mean(tempResults) | |
| } | |
| } | |
| } | |
| } | |
| #################################################################################################### | |
| # Create and test models | |
| #################################################################################################### | |
| results60<-results60[-which(results60$ACI_Deviation==0),] | |
| results12<-results12[-which(results12$ACI_Deviation==0),] | |
| attach(results12) | |
| # 12-minute ACI | |
| # Compute non-transformed global model | |
| baseModelACI12<-lmer(ACI_Deviation~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(baseModelACI12) | |
| # Check normality of residuals | |
| qqmath(baseModelACI12) | |
| shapiro.test(resid(baseModelACI12)) | |
| # Compute log-transformed global model | |
| logModelACI12<-lmer(log(ACI_Deviation)~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(logModelACI12) | |
| # Check normality of residuals | |
| qqmath(logModelACI12) | |
| shapiro.test(resid(logModelACI12)) | |
| # Compare models | |
| dredge(logModelACI12) | |
| # Create raw model with non-orthogonal polynomials and no scaling (but maintaining the log transform) | |
| logModelACI12raw<-lmer(log(ACI_Deviation)~poly(Total_Analysis_Time,3,raw=TRUE)[,1]+poly(Total_Analysis_Time,3,raw=TRUE)[,2]+poly(Total_Analysis_Time,3,raw=TRUE)[,3]+Number_of_Subdivisions+poly(Total_Analysis_Time,3,raw=TRUE)[,1]:Number_of_Subdivisions+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # 12-minute BI | |
| # Compute non-transformed global model | |
| baseModelBI12<-lmer(BI_Deviation~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(baseModelBI12) | |
| # Check normality of residuals | |
| qqmath(baseModelBI12) | |
| shapiro.test(resid(baseModelBI12)) | |
| # Compute log-transformed global model | |
| logModelBI12<-lmer(log(BI_Deviation)~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(logModelBI12) | |
| # Check normality of residuals | |
| qqmath(logModelBI12) | |
| shapiro.test(resid(logModelBI12)) | |
| # Compare models | |
| dredge(logModelBI12) | |
| # Create raw model with non-orthogonal polynomials and no scaling (but maintaining the log transform) | |
| logModelBI12raw<-lmer(log(BI_Deviation)~poly(Total_Analysis_Time,3,raw=TRUE)[,1]+poly(Total_Analysis_Time,3,raw=TRUE)[,2]+poly(Total_Analysis_Time,3,raw=TRUE)[,3]+Number_of_Subdivisions+poly(Total_Analysis_Time,3,raw=TRUE)[,1]:Number_of_Subdivisions+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # 12-minute Hf | |
| # Compute non-transformed global model | |
| baseModelHf12<-lmer(Hf_Deviation~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(baseModelHf12) | |
| # Check normality of residuals | |
| qqmath(baseModelHf12) | |
| shapiro.test(resid(baseModelHf12)) | |
| # Compute log-transformed global model | |
| logModelHf12<-lmer(log(Hf_Deviation)~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(logModelHf12) | |
| # Check normality of residuals | |
| qqmath(logModelHf12) | |
| shapiro.test(resid(logModelHf12)) | |
| # Compare models | |
| dredge(logModelHf12) | |
| # Create raw model with non-orthogonal polynomials and no scaling (but maintaining the log transform) | |
| logModelHf12raw<-lmer(log(Hf_Deviation)~poly(Total_Analysis_Time,3,raw=TRUE)[,1]+poly(Total_Analysis_Time,3,raw=TRUE)[,2]+poly(Total_Analysis_Time,3,raw=TRUE)[,3]+Number_of_Subdivisions+poly(Total_Analysis_Time,3,raw=TRUE)[,1]:Number_of_Subdivisions+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| detach(results12) | |
| attach(results60) | |
| # 60-minute ACI | |
| # Compute non-transformed global model | |
| baseModelACI60<-lmer(ACI_Deviation~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(baseModelACI60) | |
| # Check normality of residuals | |
| qqmath(baseModelACI60) | |
| shapiro.test(resid(baseModelACI60)) | |
| # Compute log-transformed global model | |
| logModelACI60<-lmer(log(ACI_Deviation)~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(logModelACI60) | |
| # Check normality of residuals | |
| qqmath(logModelACI60) | |
| shapiro.test(resid(logModelACI60)) | |
| # Compare models | |
| dredge(logModelACI60) | |
| # Create raw model with non-orthogonal polynomials and no scaling (but maintaining the log transform) | |
| logModelACI60raw<-lmer(log(ACI_Deviation)~poly(Total_Analysis_Time,3,raw=TRUE)[,1]+poly(Total_Analysis_Time,3,raw=TRUE)[,2]+poly(Total_Analysis_Time,3,raw=TRUE)[,3]+Number_of_Subdivisions+poly(Total_Analysis_Time,3,raw=TRUE)[,1]:Number_of_Subdivisions+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # 60-minute BI | |
| # Compute non-transformed global model | |
| baseModelBI60<-lmer(BI_Deviation~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(baseModelBI60) | |
| # Check normality of residuals | |
| qqmath(baseModelBI60) | |
| shapiro.test(resid(baseModelBI60)) | |
| # Compute log-transformed global model | |
| logModelBI60<-lmer(log(BI_Deviation)~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(logModelBI60) | |
| # Check normality of residuals | |
| qqmath(logModelBI60) | |
| shapiro.test(resid(logModelBI60)) | |
| # Compare models | |
| dredge(logModelBI60) | |
| # Create raw model with non-orthogonal polynomials and no scaling (but maintaining the log transform) | |
| logModelBI60raw<-lmer(log(BI_Deviation)~poly(Total_Analysis_Time,3,raw=TRUE)[,1]+poly(Total_Analysis_Time,3,raw=TRUE)[,2]+poly(Total_Analysis_Time,3,raw=TRUE)[,3]+Number_of_Subdivisions+poly(Total_Analysis_Time,3,raw=TRUE)[,1]:Number_of_Subdivisions+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # 60-minute Hf | |
| # Compute non-transformed global model | |
| baseModelHf60<-lmer(Hf_Deviation~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(baseModelHf60) | |
| # Check normality of residuals | |
| qqmath(baseModelHf60) | |
| shapiro.test(resid(baseModelHf60)) | |
| # Compute log-transformed global model | |
| logModelHf60<-lmer(log(Hf_Deviation)~scale(poly(Total_Analysis_Time,3)[,1])+scale(poly(Total_Analysis_Time,3)[,2])+scale(poly(Total_Analysis_Time,3)[,3])+scale(Number_of_Subdivisions)+scale(poly(Total_Analysis_Time,3)[,1]):scale(Number_of_Subdivisions)+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| # Check homogeneity of variance | |
| plot(logModelHf60) | |
| # Check normality of residuals | |
| qqmath(logModelHf60) | |
| shapiro.test(resid(logModelHf60)) | |
| # Compare models | |
| dredge(logModelHf60) | |
| # Create raw model with non-orthogonal polynomials and no scaling (but maintaining the log transform) | |
| logModelHf60raw<-lmer(log(Hf_Deviation)~poly(Total_Analysis_Time,3,raw=TRUE)[,1]+poly(Total_Analysis_Time,3,raw=TRUE)[,2]+poly(Total_Analysis_Time,3,raw=TRUE)[,3]+Number_of_Subdivisions+poly(Total_Analysis_Time,3,raw=TRUE)[,1]:Number_of_Subdivisions+Time_of_Day+(1|Location),REML=FALSE,na.action=na.fail) | |
| detach(results60) | |
| #################################################################################################### | |
| # Create graphics | |
| #################################################################################################### | |
| rawModels<-array(data=list(coef(logModelACI12raw)$Location,coef(logModelBI12raw)$Location,coef(logModelHf12raw)$Location,coef(logModelACI60raw)$Location,coef(logModelBI60raw)$Location,coef(logModelHf60raw)$Location,"ACI","BI","Hf","ACI","BI","Hf",12,12,12,60,60,60),dim=c(6,3),dimnames=list(NULL,c("Model","Index","Duration"))) | |
| for(i in 1:6) | |
| { | |
| # Obtain necessary values for plotting | |
| rawCoefficients<-rawModels[i,1]$Model | |
| index<-rawModels[i,2]$Index | |
| duration<-rawModels[i,3]$Duration | |
| dawnCoefficients<-rawCoefficients[,1] | |
| daytimeCoefficients<-rawCoefficients[,1]+rawCoefficients[1,6] | |
| duskCoefficients<-rawCoefficients[,1]+rawCoefficients[1,7] | |
| nighttimeCoefficients<-rawCoefficients[,1]+rawCoefficients[1,8] | |
| totalAnalysisTime1Coefficient<-rawCoefficients[1,2] | |
| totalAnalysisTime2Coefficient<-rawCoefficients[1,3] | |
| totalAnalysisTime3Coefficient<-rawCoefficients[1,4] | |
| numberOfSubdivisionsCoefficient<-rawCoefficients[1,5] | |
| interactionCoefficient<-rawCoefficients[1,9] | |
| if(duration==12) totalAnalysisTimes<-seq(0,12,1) | |
| if(duration==60) totalAnalysisTimes<-seq(0,60,1) | |
| numbersOfSubdivisions<-factors(duration)[-length(factors(duration))] | |
| dawnPredictedValues<-array(data=NA,dim=c(8,duration+1,length(factors(duration))-1),dimnames=list(c("site1","site2","site3","site4","site5","site6","site7","site8"),totalAnalysisTimes,factors(duration)[-length(factors(duration))])) | |
| daytimePredictedValues<-dawnPredictedValues | |
| duskPredictedValues<-dawnPredictedValues | |
| nighttimePredictedValues<-dawnPredictedValues | |
| meanPredictedValuesAllSites<-dawnPredictedValues | |
| for(j in 1:8) | |
| { | |
| for(k in 1:(length(factors(duration))-1)) | |
| { | |
| numberOfSubdivisions<-numbersOfSubdivisions[k] | |
| dawnPredictedValues[j,,k]<-exp(dawnCoefficients[j]+totalAnalysisTime1Coefficient*totalAnalysisTimes+totalAnalysisTime2Coefficient*(totalAnalysisTimes^2)+totalAnalysisTime3Coefficient*(totalAnalysisTimes^3)+numberOfSubdivisionsCoefficient*numberOfSubdivisions+interactionCoefficient*totalAnalysisTimes*numberOfSubdivisions) | |
| daytimePredictedValues[j,,k]<-exp(daytimeCoefficients[j]+totalAnalysisTime1Coefficient*totalAnalysisTimes+totalAnalysisTime2Coefficient*(totalAnalysisTimes^2)+totalAnalysisTime3Coefficient*(totalAnalysisTimes^3)+numberOfSubdivisionsCoefficient*numberOfSubdivisions+interactionCoefficient*totalAnalysisTimes*numberOfSubdivisions) | |
| duskPredictedValues[j,,k]<-exp(duskCoefficients[j]+totalAnalysisTime1Coefficient*totalAnalysisTimes+totalAnalysisTime2Coefficient*(totalAnalysisTimes^2)+totalAnalysisTime3Coefficient*(totalAnalysisTimes^3)+numberOfSubdivisionsCoefficient*numberOfSubdivisions+interactionCoefficient*totalAnalysisTimes*numberOfSubdivisions) | |
| nighttimePredictedValues[j,,k]<-exp(nighttimeCoefficients[j]+totalAnalysisTime1Coefficient*totalAnalysisTimes+totalAnalysisTime2Coefficient*(totalAnalysisTimes^2)+totalAnalysisTime3Coefficient*(totalAnalysisTimes^3)+numberOfSubdivisionsCoefficient*numberOfSubdivisions+interactionCoefficient*totalAnalysisTimes*numberOfSubdivisions) | |
| } | |
| } | |
| meanDawnPredictedValues<-t(colMeans(dawnPredictedValues)) | |
| meanDaytimePredictedValues<-t(colMeans(daytimePredictedValues)) | |
| meanDuskPredictedValues<-t(colMeans(duskPredictedValues)) | |
| meanNighttimePredictedValues<-t(colMeans(nighttimePredictedValues)) | |
| meanPredictedValuesAllSites<-(dawnPredictedValues+daytimePredictedValues+duskPredictedValues+nighttimePredictedValues)/4 | |
| meanPredictedValues<-t(colMeans(meanPredictedValuesAllSites)) | |
| upper95Values<-matrix(NA,(length(numbersOfSubdivisions)-1),(duration+1)) | |
| lower95Values<-upper95Values | |
| for(j in 1:(length(numbersOfSubdivisions)-1)) upper95Values[j,]<-meanPredictedValues[j,]+1.96*apply(meanPredictedValuesAllSites[,,j],2,sd)/sqrt(8) | |
| for(j in 1:(length(numbersOfSubdivisions)-1)) lower95Values[j,]<-meanPredictedValues[j,]-1.96*apply(meanPredictedValuesAllSites[,,j],2,sd)/sqrt(8) | |
| colorSequence<-c("#008000","#d95f02","#7570b3","#e7298a","#7dc17d","#e6ab02","#a6761d","#666666") | |
| if(index=="ACI"|index=="BI") ymax<-trunc(max(upper95Values)+1) | |
| if(index=="Hf") ymax<-round(max(upper95Values)+0.005,2) | |
| # Plot 1-subdivision values for times of day | |
| pdf(paste("/Users/dantefrancomano/Documents/Purdue/Projects/Temporal_Variability_Duty_Cycle/Figures/Deviation/",index,"_",duration,"_Times_of_Day.pdf",sep=""),5,4,pointsize=8) | |
| plot(1:2,1:2,type="n",xlim=c(0,duration),ylim=c(0,ymax),axes=FALSE,main=paste(index," ",duration,"-minute Window",sep=""),col.main="grey20",xlab="",ylab="") | |
| if(index=="ACI"|index=="BI") abline(h=seq(0,ymax,ymax/5),col="grey90") | |
| if(index=="Hf") abline(h=seq(0,ymax,ymax/6),col="grey90") | |
| abline(v=seq(0,duration,duration/6),col="grey90") | |
| box(col="grey30") | |
| axis(1,col="grey30",col.axis="grey30") | |
| axis(2,col="grey30",col.axis="grey30") | |
| mtext("Total Analysis Time (minutes)",side=1,col="grey30",padj=4) | |
| mtext(paste(index,"Deviation"),side=2,col="grey30",padj=-4) | |
| lines(totalAnalysisTimes,meanDawnPredictedValues[1,],col=colorSequence[1],lwd=3) | |
| lines(totalAnalysisTimes,meanDaytimePredictedValues[1,],col=colorSequence[2],lwd=3) | |
| lines(totalAnalysisTimes,meanDuskPredictedValues[1,],col=colorSequence[3],lwd=3) | |
| lines(totalAnalysisTimes,meanNighttimePredictedValues[1,],col=colorSequence[6],lwd=3) | |
| legend("topright",inset=0.025,legend=c("Dawn","Daytime","Dusk","Nighttime"),lty=1,lwd=3,col=colorSequence[c(1,2,3,6)],text.col="grey30",bg="white",box.col="grey30") | |
| dev.off() | |
| # Plot 1-subdivision global mean and confidence intervals | |
| pdf(paste("/Users/dantefrancomano/Documents/Purdue/Projects/Temporal_Variability_Duty_Cycle/Figures/Deviation/",index,"_",duration,"_Confidence_Interval.pdf",sep=""),5,4,pointsize=8) | |
| plot(1:2,1:2,type="n",xlim=c(0,duration),ylim=c(0,ymax),axes=FALSE,main=paste(index," ",duration,"-minute Window",sep=""),col.main="grey20",xlab="",ylab="") | |
| if(index=="ACI"|index=="BI") abline(h=seq(0,ymax,ymax/5),col="grey90") | |
| if(index=="Hf") abline(h=seq(0,ymax,ymax/6),col="grey90") | |
| abline(v=seq(0,duration,duration/6),col="grey90") | |
| box(col="grey30") | |
| axis(1,col="grey30",col.axis="grey30") | |
| axis(2,col="grey30",col.axis="grey30") | |
| mtext("Total Analysis Time (minutes)",side=1,col="grey30",padj=4) | |
| mtext(paste(index,"Deviation"),side=2,col="grey30",padj=-4) | |
| polygon(c(totalAnalysisTimes,rev(totalAnalysisTimes)),c(upper95Values[1,],rev(lower95Values[1,])),density=1000,col=rgb(t(col2rgb("grey30")/255),alpha=0.05)) | |
| lines(totalAnalysisTimes,meanPredictedValues[1,],lwd=3,col="grey30") | |
| legend("topright",inset=0.025,legend=c("Mean","95% Confidence Interval"),lty=c(1,0),lwd=c(3,0),pch=c(NA,15),pt.cex=2.5,col=c("grey30",rgb(t(col2rgb("grey30")/255),alpha=0.5)),text.col="grey30",bg="white",box.col="grey30") | |
| dev.off() | |
| # Plot multiple subdivisions | |
| pdf(paste("/Users/dantefrancomano/Documents/Purdue/Projects/Temporal_Variability_Duty_Cycle/Figures/Deviation/",index,"_",duration,"_Confidence_Interval_Subdivisions_Lines.pdf",sep=""),5,4,pointsize=8) | |
| plot(1:2,1:2,type="n",xlim=c(0,duration),ylim=c(0,ymax),axes=FALSE,main=paste(index," ",duration,"-minute window",sep=""),col.main="grey20",xlab="",ylab="") | |
| if(index=="ACI"|index=="BI") abline(h=seq(0,ymax,ymax/5),col="grey90") | |
| if(index=="Hf") abline(h=seq(0,ymax,ymax/6),col="grey90") | |
| abline(v=seq(0,duration,duration/6),col="grey90") | |
| box(col="grey30") | |
| axis(1,col="grey30",col.axis="grey30") | |
| axis(2,col="grey30",col.axis="grey30") | |
| mtext("total analysis time (minutes)",side=1,col="grey30",padj=4) | |
| mtext(paste(index," deviation (absolute value)",sep=""),side=2,col="grey30",padj=-4) | |
| polygon(c(totalAnalysisTimes,rev(totalAnalysisTimes)),c(upper95Values[1,],rev(lower95Values[1,])),density=1000,col=rgb(t(col2rgb("grey80")/255),alpha=0.05)) | |
| if(duration==12) | |
| { | |
| lines(totalAnalysisTimes,meanPredictedValues[5,],lwd=3,lty=1,col="grey70") | |
| lines(totalAnalysisTimes,meanPredictedValues[4,],lwd=3,lty=1,col="grey60") | |
| lines(totalAnalysisTimes,meanPredictedValues[3,],lwd=3,lty=1,col="grey50") | |
| lines(totalAnalysisTimes,meanPredictedValues[2,],lwd=3,lty=1,col="grey40") | |
| lines(totalAnalysisTimes,meanPredictedValues[1,],lwd=3,lty=1,col="grey30") | |
| legend("topright",inset=0.025,legend=c("1 subdivision mean","95% confidence interval","2 subdivision mean","3 subdivision mean","4 subdivision mean","6 subdivision mean"),lty=c(1,0,1,1,1,1),lwd=3,pch=c(NA,15,NA,NA,NA,NA),pt.cex=2.5,col=c("grey30",rgb(t(col2rgb("grey80")/255),alpha=0.5),"grey40","grey50","grey60","grey70"),text.col="grey30",bg="white",box.col="grey30") | |
| }else{ | |
| lines(totalAnalysisTimes,meanPredictedValues[9,],lwd=3,lty=1,col="grey80") | |
| lines(totalAnalysisTimes,meanPredictedValues[7,],lwd=3,lty=1,col="grey70") | |
| lines(totalAnalysisTimes,meanPredictedValues[6,],lwd=3,lty=1,col="grey60") | |
| lines(totalAnalysisTimes,meanPredictedValues[4,],lwd=3,lty=1,col="grey50") | |
| lines(totalAnalysisTimes,meanPredictedValues[2,],lwd=3,lty=1,col="grey40") | |
| lines(totalAnalysisTimes,meanPredictedValues[1,],lwd=3,lty=1,col="grey30") | |
| legend("topright",inset=0.025,legend=c("1 subdivision mean","95% confidence interval","2 subdivision mean","4 subdivision mean","6 subdivision mean","10 subdivision mean","15 subdivision mean"),lty=c(1,0,1,1,1,1,1),lwd=3,pch=c(NA,15,NA,NA,NA,NA,NA),pt.cex=2.5,col=c("grey30",rgb(t(col2rgb("grey80")/255),alpha=0.5),"grey40","grey50","grey60","grey70","grey80"),text.col="grey30",bg="white",box.col="grey30") | |
| } | |
| dev.off() | |
| pdf(paste("/Users/dantefrancomano/Documents/Purdue/Projects/Temporal_Variability_Duty_Cycle/Figures/Deviation/",index,"_",duration,"_Confidence_Interval_Subdivisions_Text.pdf",sep=""),5,4,pointsize=8) | |
| plot(1:2,1:2,type="n",xlim=c(0,duration),ylim=c(0,ymax),axes=FALSE,main=paste(index," ",duration,"-minute Window",sep=""),col.main="grey20",xlab="",ylab="") | |
| if(index=="ACI"|index=="BI") abline(h=seq(0,ymax,ymax/5),col="grey90") | |
| if(index=="Hf") abline(h=seq(0,ymax,ymax/6),col="grey90") | |
| abline(v=seq(0,duration,duration/6),col="grey90") | |
| box(col="grey30") | |
| axis(1,col="grey30",col.axis="grey30") | |
| axis(2,col="grey30",col.axis="grey30") | |
| mtext("Total Analysis Time (minutes)",side=1,col="grey30",padj=4) | |
| mtext(paste(index,"Deviation"),side=2,col="grey30",padj=-4) | |
| polygon(c(totalAnalysisTimes,rev(totalAnalysisTimes)),c(upper95Values[1,],rev(lower95Values[1,])),density=1000,col=rgb(t(col2rgb("grey30")/255),alpha=0.05)) | |
| lines(totalAnalysisTimes,meanPredictedValues[1,],lwd=3,lty=1,col="grey30") | |
| if(duration==12) | |
| { | |
| text(x=c(2,4,6),y=meanPredictedValues[2,c(3,5,7)],labels="2",font=2,col="grey30") | |
| text(x=c(3,6,9),y=meanPredictedValues[3,c(4,7,10)],labels="3",font=2,col="grey30") | |
| text(x=c(4,8),y=meanPredictedValues[4,c(5,9)],labels="4",font=2,col="grey30") | |
| text(x=6,y=meanPredictedValues[5,7],labels="6",font=2,col="grey30") | |
| legend("topright",inset=0.025,legend=c("1 Subdivision Mean","95% Confidence Interval","Number of Subdivisions"),lty=c(1,0,0),lwd=3,pch=c(NA,15,NA),pt.cex=2.5,col=c("grey30",rgb(t(col2rgb("grey30")/255),alpha=0.5),"grey30"),text.col="grey30",bg="white",box.col="grey30") | |
| legend("topright",inset=c(0.028,0.025),legend=c("","","# "),text.col="grey30",text.font=2,bg="transparent",box.col=NA) | |
| }else{ | |
| text(x=seq(2,58,2),y=meanPredictedValues[2,seq(3,59,2)],labels="2",font=2,col="grey30") | |
| text(x=seq(4,56,4),y=meanPredictedValues[4,seq(5,57,4)],labels="4",font=2,col="grey30") | |
| text(x=seq(6,54,6),y=meanPredictedValues[6,seq(7,55,6)],labels="6",font=2,col="grey30") | |
| text(x=seq(12,48,12),y=meanPredictedValues[8,seq(13,49,12)],labels="12",font=2,col="grey30") | |
| legend("topright",inset=0.025,legend=c("1 Subdivision Mean","95% Confidence Interval","Number of Subdivisions"),lty=c(1,0,0),lwd=3,pch=c(NA,15,NA),pt.cex=2.5,col=c("grey30",rgb(t(col2rgb("grey30")/255),alpha=0.5),"grey30"),text.col="grey30",bg="white",box.col="grey30") | |
| legend("topright",inset=c(0.028,0.025),legend=c("","","# "),text.col="grey30",text.font=2,bg="transparent",box.col=NA) | |
| } | |
| dev.off() | |
| } | |
| #################################################################################################### | |
| # Create combined plot | |
| #################################################################################################### | |
| pdf("/Users/dantefrancomano/Documents/Purdue/Projects/Temporal_Variability_Duty_Cycle/Figures/Question_2_Combined_Plots.pdf",7.5,10,pointsize=8) | |
| par(mfcol=c(3,2)) | |
| par(oma=c(4,4,0,0)) | |
| par(mar=c(2,4,4,2)+0.1) | |
| letters<-c("A","B","C","D","E","F") | |
| for(i in 1:6) | |
| { | |
| # Obtain necessary values for plotting | |
| rawCoefficients<-rawModels[i,1]$Model | |
| index<-rawModels[i,2]$Index | |
| duration<-rawModels[i,3]$Duration | |
| dawnCoefficients<-rawCoefficients[,1] | |
| daytimeCoefficients<-rawCoefficients[,1]+rawCoefficients[1,6] | |
| duskCoefficients<-rawCoefficients[,1]+rawCoefficients[1,7] | |
| nighttimeCoefficients<-rawCoefficients[,1]+rawCoefficients[1,8] | |
| totalAnalysisTime1Coefficient<-rawCoefficients[1,2] | |
| totalAnalysisTime2Coefficient<-rawCoefficients[1,3] | |
| totalAnalysisTime3Coefficient<-rawCoefficients[1,4] | |
| numberOfSubdivisionsCoefficient<-rawCoefficients[1,5] | |
| interactionCoefficient<-rawCoefficients[1,9] | |
| if(duration==12) totalAnalysisTimes<-seq(0,12,1) | |
| if(duration==60) totalAnalysisTimes<-seq(0,60,1) | |
| numbersOfSubdivisions<-factors(duration)[-length(factors(duration))] | |
| dawnPredictedValues<-array(data=NA,dim=c(8,duration+1,length(factors(duration))-1),dimnames=list(c("site1","site2","site3","site4","site5","site6","site7","site8"),totalAnalysisTimes,factors(duration)[-length(factors(duration))])) | |
| daytimePredictedValues<-dawnPredictedValues | |
| duskPredictedValues<-dawnPredictedValues | |
| nighttimePredictedValues<-dawnPredictedValues | |
| meanPredictedValuesAllSites<-dawnPredictedValues | |
| for(j in 1:8) | |
| { | |
| for(k in 1:(length(factors(duration))-1)) | |
| { | |
| numberOfSubdivisions<-numbersOfSubdivisions[k] | |
| dawnPredictedValues[j,,k]<-exp(dawnCoefficients[j]+totalAnalysisTime1Coefficient*totalAnalysisTimes+totalAnalysisTime2Coefficient*(totalAnalysisTimes^2)+totalAnalysisTime3Coefficient*(totalAnalysisTimes^3)+numberOfSubdivisionsCoefficient*numberOfSubdivisions+interactionCoefficient*totalAnalysisTimes*numberOfSubdivisions) | |
| daytimePredictedValues[j,,k]<-exp(daytimeCoefficients[j]+totalAnalysisTime1Coefficient*totalAnalysisTimes+totalAnalysisTime2Coefficient*(totalAnalysisTimes^2)+totalAnalysisTime3Coefficient*(totalAnalysisTimes^3)+numberOfSubdivisionsCoefficient*numberOfSubdivisions+interactionCoefficient*totalAnalysisTimes*numberOfSubdivisions) | |
| duskPredictedValues[j,,k]<-exp(duskCoefficients[j]+totalAnalysisTime1Coefficient*totalAnalysisTimes+totalAnalysisTime2Coefficient*(totalAnalysisTimes^2)+totalAnalysisTime3Coefficient*(totalAnalysisTimes^3)+numberOfSubdivisionsCoefficient*numberOfSubdivisions+interactionCoefficient*totalAnalysisTimes*numberOfSubdivisions) | |
| nighttimePredictedValues[j,,k]<-exp(nighttimeCoefficients[j]+totalAnalysisTime1Coefficient*totalAnalysisTimes+totalAnalysisTime2Coefficient*(totalAnalysisTimes^2)+totalAnalysisTime3Coefficient*(totalAnalysisTimes^3)+numberOfSubdivisionsCoefficient*numberOfSubdivisions+interactionCoefficient*totalAnalysisTimes*numberOfSubdivisions) | |
| } | |
| } | |
| meanDawnPredictedValues<-t(colMeans(dawnPredictedValues)) | |
| meanDaytimePredictedValues<-t(colMeans(daytimePredictedValues)) | |
| meanDuskPredictedValues<-t(colMeans(duskPredictedValues)) | |
| meanNighttimePredictedValues<-t(colMeans(nighttimePredictedValues)) | |
| meanPredictedValuesAllSites<-(dawnPredictedValues+daytimePredictedValues+duskPredictedValues+nighttimePredictedValues)/4 | |
| meanPredictedValues<-t(colMeans(meanPredictedValuesAllSites)) | |
| upper95Values<-matrix(NA,(length(numbersOfSubdivisions)-1),(duration+1)) | |
| lower95Values<-upper95Values | |
| for(j in 1:(length(numbersOfSubdivisions)-1)) upper95Values[j,]<-meanPredictedValues[j,]+1.96*apply(meanPredictedValuesAllSites[,,j],2,sd)/sqrt(8) | |
| for(j in 1:(length(numbersOfSubdivisions)-1)) lower95Values[j,]<-meanPredictedValues[j,]-1.96*apply(meanPredictedValuesAllSites[,,j],2,sd)/sqrt(8) | |
| colorSequence<-c("#008000","#d95f02","#7570b3","#e7298a","#7dc17d","#e6ab02","#a6761d","#666666") | |
| if(index=="ACI"|index=="BI") ymax<-trunc(max(upper95Values)+1) | |
| if(index=="Hf") ymax<-round(max(upper95Values)+0.005,2) | |
| plot(1:2,1:2,type="n",xlim=c(0,duration),ylim=c(0,ymax),axes=FALSE,xlab="",ylab="") | |
| if(index=="ACI"|index=="BI") | |
| { | |
| if(index=="ACI") abline(h=seq(0,ymax,ymax/5),col="grey90") | |
| if(index=="BI") abline(h=seq(0,ymax,ymax/6),col="grey90") | |
| title(paste(letters[i],". ",index,": ",duration,"-minute full file length",sep=""),adj=0,cex.main=2,col.main="grey20") | |
| }else{ | |
| abline(h=seq(0,ymax,ymax/6),col="grey90") | |
| if(duration==12) | |
| { | |
| title(expression(bold("C. H"["f"]*": 12-minute full file length")),adj=0,cex.main=2,col.main="grey20") | |
| }else{ | |
| title(expression(bold("F. H"["f"]*": 60-minute full file length")),adj=0,cex.main=2,col.main="grey20") | |
| } | |
| } | |
| abline(v=seq(0,duration,duration/6),col="grey90") | |
| box(col="grey30") | |
| axis(1,col="grey30",col.axis="grey30",cex.axis=2,padj=0.5) | |
| axis(2,col="grey30",col.axis="grey30",cex.axis=2,las=2) | |
| polygon(c(totalAnalysisTimes,rev(totalAnalysisTimes)),c(upper95Values[1,],rev(lower95Values[1,])),density=1000,col=rgb(t(col2rgb("grey80")/255),alpha=0.05)) | |
| if(duration==12) | |
| { | |
| lines(totalAnalysisTimes,meanPredictedValues[5,],lwd=3,lty=1,col="grey70") | |
| lines(totalAnalysisTimes,meanPredictedValues[4,],lwd=3,lty=1,col="grey60") | |
| lines(totalAnalysisTimes,meanPredictedValues[3,],lwd=3,lty=1,col="grey50") | |
| lines(totalAnalysisTimes,meanPredictedValues[2,],lwd=3,lty=1,col="grey40") | |
| lines(totalAnalysisTimes,meanPredictedValues[1,],lwd=3,lty=1,col="grey30") | |
| legend("topright",inset=0.025,legend=c("1 subdivision mean predicted values","1 subdivision 95% confidence interval","2 subdivision mean predicted values","3 subdivision mean predicted values","4 subdivision mean predicted values","6 subdivision mean predicted values"),lty=c(1,0,1,1,1,1),lwd=3,pch=c(NA,15,NA,NA,NA,NA),pt.cex=4,col=c("grey30",rgb(t(col2rgb("grey80")/255),alpha=0.5),"grey40","grey50","grey60","grey70"),text.col="grey30",bg="white",box.col="grey30",cex=1.5) | |
| }else{ | |
| lines(totalAnalysisTimes,meanPredictedValues[8,],lwd=3,lty=1,col="grey70") | |
| lines(totalAnalysisTimes,meanPredictedValues[6,],lwd=3,lty=1,col="grey60") | |
| lines(totalAnalysisTimes,meanPredictedValues[4,],lwd=3,lty=1,col="grey50") | |
| lines(totalAnalysisTimes,meanPredictedValues[2,],lwd=3,lty=1,col="grey40") | |
| lines(totalAnalysisTimes,meanPredictedValues[1,],lwd=3,lty=1,col="grey30") | |
| legend("topright",inset=0.025,legend=c("1 subdivision mean predicted values","1 subdivision 95% confidence interval","2 subdivision mean predicted values","4 subdivision mean predicted values","6 subdivision mean predicted values","12 subdivision mean predicted values"),lty=c(1,0,1,1,1,1),lwd=3,pch=c(NA,15,NA,NA,NA,NA),pt.cex=4,col=c("grey30",rgb(t(col2rgb("grey80")/255),alpha=0.5),"grey40","grey50","grey60","grey70"),text.col="grey30",bg="white",box.col="grey30",cex=1.5) | |
| } | |
| } | |
| mtext("total analysis time (minutes)",side=1,outer=TRUE,col="grey30",cex=2,padj=1.4) | |
| mtext("subsample error (absolute value)",side=2,outer=TRUE,col="grey30",cex=2,padj=-0.3) | |
| dev.off() | |
| par(mfrow=c(1,1)) | |
| par(oma=c(1,1,1,1)) | |
| par(mar=c(5,4,4,2)+0.1) |