library("ggplot2") library("glmmADMB") library("plyr") library("reshape2") #load saved simulation data load("srs_glmm_results_total_asco_meanCenter.rdata") load("MMNND_glmm_results_total_asco_meanCenter.rdata") load("GLMM_total_asco_sim.rdata") params_MMNND_total_asco<-matrix(nrow=1000,ncol=6) for (i in 1:length(MMNND_glmm_results_total_asco_meanCenter)){ for (j in 1:6){ MMNND_glmm_results_total_asco_meanCenter[[i]]$b[[j]] -> params_MMNND_total_asco[i,j] #print(i) } } params_MMNND_total_asco_dframe <- as.data.frame(params_MMNND_total_asco) colnames(params_MMNND_total_asco_dframe)<-c("Intercept", "geo_dis", "gen_dis", "geo_sq", "gen_dis_sq", "geo:gen_dis") params_MMNND_total_asco_dframe<-melt(params_MMNND_total_asco_dframe) params_MMNND_total_asco_dframe$sampling_method <- "SPREAD" #for SRS params_srs_total_asco<-matrix(nrow=1000,ncol=6) for (i in 1:length(srs_glmm_results_total_asco_meanCenter)){ for (j in 1:6){ srs_glmm_results_total_asco_meanCenter[[i]]$b[[j]] -> params_srs_total_asco[i,j] #print(i) } } # for srs params_srs_total_asco_dframe <- as.data.frame(params_srs_total_asco) colnames(params_srs_total_asco_dframe)<-c("Intercept", "geo_dis", "gen_dis", "geo_sq", "gen_dis_sq", "geo:gen_dis") params_srs_total_asco_dframe<-melt(params_srs_total_asco_dframe) params_srs_total_asco_dframe$sampling_method <- "SRS" SRS_MMNND_dframe_params <- rbind(params_MMNND_total_asco_dframe, params_srs_total_asco_dframe) SRS_MMNND_dframe_params$variable<- as.factor(SRS_MMNND_dframe_params$variable) colnames(SRS_MMNND_dframe_params)<-c("Parameter", "Parameter_Value", "Sampling_Method") actual_params<-data.frame(Parameter=as.factor(c("Intercept", "geo_dis", "gen_dis", "geo_sq", "gen_dis_sq", "geo:gen_dis")), Parameter_Value_original=c(MP1_total_asco_sim$b[[1]], MP1_total_asco_sim$b[[2]], MP1_total_asco_sim$b[[3]], MP1_total_asco_sim$b[[4]], MP1_total_asco_sim$b[[5]], MP1_total_asco_sim$b[[6]]), Std_error=c(stdEr(MP1_total_asco_sim)[[1]], stdEr(MP1_total_asco_sim)[[2]], stdEr(MP1_total_asco_sim)[[3]], stdEr(MP1_total_asco_sim)[[4]], stdEr(MP1_total_asco_sim)[[5]], stdEr(MP1_total_asco_sim)[[6]])) #violin plots of parameter values (Figure 3A) params_violin_plot <- ggplot(SRS_MMNND_dframe_params, aes(x=Sampling_Method, y=Parameter_Value)) + geom_violin() + facet_wrap(~Parameter, scales="free_y", ncol=6) + geom_hline(data = actual_params, aes(yintercept = Parameter_Value_original, colour="black"))+ theme_bw() + theme( strip.background=element_rect(colour = NA, fill = NA))+ xlab("Sampling Method") + ylab("Parameter Value") ggsave(filename = "Figure_3A.pdf", plot = params_violin_plot, height=5, width=10) stderror_MMNND_total_asco<-matrix(nrow=1000,ncol=6) for (i in 1:length(MMNND_glmm_results_total_asco_meanCenter)){ for (j in 1:6){ stdEr(MMNND_glmm_results_total_asco_meanCenter[[i]])[[j]] -> stderror_MMNND_total_asco[i,j] } } # for MMNND stderror_MMNND_total_asco_dframe <- as.data.frame(stderror_MMNND_total_asco) colnames(stderror_MMNND_total_asco_dframe)<-c("Intercept", "geo_dis", "gen_dis", "geo_sq", "gen_dis_sq", "geo:gen_dis") stderror_MMNND_total_asco_dframe<-melt(stderror_MMNND_total_asco_dframe) stderror_MMNND_total_asco_dframe$sampling_method <- "SPREAD" #for SRS stderror_srs_total_asco<-matrix(nrow=1000,ncol=6) for (i in 1:length(srs_glmm_results_total_asco_meanCenter)){ for (j in 1:6){ stdEr(srs_glmm_results_total_asco_meanCenter[[i]])[[j]] -> stderror_srs_total_asco[i,j] #print(i) } } # for srs stderror_srs_total_asco_dframe <- as.data.frame(stderror_srs_total_asco) colnames(stderror_srs_total_asco_dframe)<-c("Intercept", "geo_dis", "gen_dis", "geo_sq", "gen_dis_sq", "geo:gen_dis") stderror_srs_total_asco_dframe<-melt(stderror_srs_total_asco_dframe) stderror_srs_total_asco_dframe$sampling_method <- "SRS" #combine data frames for comparative analysis: SRS_MMNND_dframe_std_error <- rbind(stderror_MMNND_total_asco_dframe, stderror_srs_total_asco_dframe) SRS_MMNND_dframe_std_error$variable<- as.factor(SRS_MMNND_dframe_std_error$variable) colnames(SRS_MMNND_dframe_std_error)<-c("Parameter", "Standard_Error", "Sampling_Method") # Dataframe of standard errors (Table 2) summary_std_errors<-ddply(SRS_MMNND_dframe_std_error, .(Sampling_Method, Parameter), summarize, mean_std_error=mean(Standard_Error), variance_std_error = var(Standard_Error)) write.csv(summary_std_errors, "table2.csv") #violin plots of standard errors (Figure 3B) stderror_violin_plot <- ggplot(SRS_MMNND_dframe_std_error, aes(x=Sampling_Method, y=Standard_Error)) + geom_violin() + facet_wrap(~Parameter, scales="free_y", ncol=6) + theme_bw() + theme(legend.position = "none", strip.background=element_rect(colour = NA, fill = NA))+ xlab("Sampling Method") + ylab("Standard Error") #save plot ggsave(filename = "Figure_3B.pdf", plot = stderror_violin_plot, height=5, width=10)