.Rhistory

ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred))), color="black") +
labs(title=paste("Population Abundance Prediction"), x="X (30m increment)", y="Y (30m increment)", fill="Cover")
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance"), x="X (30m increment)", y="Y (30m increment)", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA)
ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred))), color="black") +
labs(title=paste("Population Abundance Prediction"), x="X (30m increment)", y="Y (30m increment)", fill="Cover")
min(test$easting)
max(test$easting)
min(test$northing)
max(test$northing)
4303800
ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred))), color="black") +
labs(title=paste("Population Abundance Prediction"), x="X (30m increment)", y="Y (30m increment)", fill="Cover") +
scale_x_continuous(limits = c(365430, 366240)) +
scale_y_continuous(limits = c(4303800, 4304430))
ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred))), color="black") +
labs(title=paste("Population Abundance Prediction"), x="X (30m increment)", y="Y (30m increment)", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance"), x="X (30m increment)", y="Y (30m increment)", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) + +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance"), x="X (30m increment)", y="Y (30m increment)", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
# edit scale for validation
species.map[which(is.na(species.map[,1])),1] <- "none"
species.map[which(species.map[,1] == 0),1] <- "few"
species.map[which(species.map[,1] == 1),1] <- "few more"
species.map[which(species.map[,1] == 2),1] <- "little"
species.map[which(species.map[,1] == 3),1] <- "some"
species.map[which(species.map[,1] == 4),1] <- "most"
species.map[which(species.map[,1] == 5),1] <- "all"
species.map[,z] <- as.factor(species.map[,1])
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance"), x="X (30m increment)", y="Y (30m increment)", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance"), x="X (30m increment)", y="Y (30m increment)", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_discrete(limits=c("none", "few", "few more", "little", "some", "most", "all")) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
species.map$scam <- factor(species.map$scam, levels = c("none", "few", "few more", "little", "some", "most", "all"))
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance"), x="X (30m increment)", y="Y (30m increment)", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
test$pred <- factor(test$pred, levels = c("none", "few", "few more", "little", "some", "most", "all"))
test$pred <- factor(test$pred, levels = c("none", "few", "few more", "little", "some", "most"))
pred <- factor(pred, levels = c("none", "few", "few more", "little", "some", "most"))
pred <- factor(pred, levels = c("none", "few", "few more", "little", "some", "most", "all"))
ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred))), color="black") +
labs(title=paste("Population Abundance Prediction"), x="X (30m increment)", y="Y (30m increment)", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance of SERC under Landsat Plots"), x="Easting", y="Northing", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance of SERC"), x="Easting", y="Northing", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance of SERC Plots"), x="Easting", y="Northing", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred))), color="black") +
labs(title=paste("Population Abundance Prediction"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred))), color="black") +
labs(title=paste("Population Abundance Prediction"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
setwd("~/Documents/Junior_Year/DISC_REU/DISC_chesapeake/")
ggsave("plots/prediction_visual.png")
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance of SERC Plots"), x="Easting", y="Northing", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
ggsave("plots/prediction_validation.png")
# model with plot proportion as a value
model2 = randomForest(scam ~ .,data = train[,-c(1:3)], keep.forest=TRUE)
model2
plot(model2)
varImpPlot(model2)
pred2 <- predict(model2, newdata = test)
table(pred2, test$scam)
pred2 <- factor(pred, levels = c("none", "few", "few more", "little", "some", "most", "all"))
pred2.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste("Population Abundance Prediction w/ Overlap"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred2.graph
pred.graph
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance of SERC Plots"), x="Easting", y="Northing", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
model2 = randomForest(scam ~ .,data = train[,-c(1:4, 9, 10, 18)], keep.forest=TRUE)
model2
plot(model2)
varImpPlot(model2)
pred2 <- predict(model2, newdata = test)
table(pred2, test$scam)
pred2 <- factor(pred, levels = c("none", "few", "few more", "little", "some", "most", "all"))
pred2.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste("Population Abundance Prediction w/ Overlap"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred2.graph
model = randomForest(scam ~ .,data = train[,-c(1:4)], keep.forest=TRUE)
model
plot(model)
varImpPlot(model)
pred <- predict(model, newdata = test)
table(pred, test$scam)
pred <- factor(pred, levels = c("none", "few", "few more", "little", "some", "most", "all"))
pred.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred))), color="black") +
labs(title=paste("Population Abundance Prediction"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred.graph
pred2.graph
ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance of SERC Plots"), x="Easting", y="Northing", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred2.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste("Population Abundance Prediction w/ less vars"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred2.graph
pred2.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred2.graph
pred.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred))), color="black") +
labs(title=paste(curr.species, "Population Abundance Prediction"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred.graph
pred.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred))), color="black") +
labs(title=paste(curr.species, "Population Abundance Prediction w/ all vars"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred.graph
pred2.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred2.graph
ggsave("plots/prediction_visual2.png")
pred.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred))), color="black") +
labs(title=paste(curr.species, "Population Abundance Prediction w/ all vars"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred.graph
ggsave("plots/prediction_visual1.png")
par(mfrow=c(1,2))
pred.graph
pred2.graph
grid.arrange(pred.graph, pred2.graph, ncol=2)
library(gridExtra)
grid.arrange(pred.graph, pred2.graph, ncol=2)
serc.plots <- ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance of SERC Plots"), x="Easting", y="Northing", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
serc.plots
grid.arrange(serc.plots, pred2.graph, ncol=2)
ggsave("plots/plot_comparison.png")
ggsave("plots/plot_comparison.png", arrangeGrob(serc.plots, pred2.graph, ncol=2))
ggsave("plots/plot_comparison.png", arrangeGrob(serc.plots, pred2.graph, ncol=2), width = 10)
ggsave("plots/plot_comparison.png", arrangeGrob(serc.plots, pred2.graph, ncol=2), width = 11)
pred2 <- factor(pred, levels = c("none", "few", "few more", "little", "some", "most", "all"))
cols = c("none"="CCCCCC", "few" = "FFFFCC", "few more" = "CCFF99", "little" = "99FF00", "some" = "33CC33", "most" = "006600", "all" = "003300")
pred2.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470)) +
scale_colour_manual(values = cols)
pred2.graph
pred2.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470)) +
scale_colour_discrete(guide=guide_legend(override.aes=aes(values = cols)))
pred2.graph <- ggplot() +
geom_rect(data=test, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
scale_colour_discrete(guide=guide_legend(override.aes=aes(values = cols))) +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred2.graph <- ggplot(data=test) +
geom_rect(aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470))
pred2.graph <- ggplot(data=test) +
geom_rect(aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470)) +
scale_fill_manual(values = cols)
pred2.graph
pred2 <- factor(pred, levels = c("none", "few", "few more", "little", "some", "most", "all"))
cols = c("none"="red", "few" = "FFFFCC", "few more" = "CCFF99", "little" = "99FF00", "some" = "33CC33", "most" = "006600", "all" = "003300")
pred2.graph <- ggplot(data=test) +
geom_rect(aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470)) +
scale_fill_manual(values = cols)
pred2.graph
ggsave("plots/prediction_visual2.png")
pred2 <- factor(pred, levels = c("none", "few", "few more", "little", "some", "most", "all"))
cols = c("none"="#ceb467", "few" = "#ace5b2", "few more" = "#7aef87", "little" = "#57d165", "some" = "#25c637", "most" = "#197f24", "all" = "#115118")
pred2.graph <- ggplot(data=test) +
geom_rect(aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470)) +
scale_fill_manual(values = cols)
pred2.graph
pred2 <- factor(pred, levels = c("none", "few", "few more", "little", "some", "most", "all"))
cols = c("none"="#ceb467", "few" = "#ace5b2", "few more" = "#7aef87", "little" = "#57d165", "some" = "#1f961b", "most" = "#197f24", "all" = "#115118")
pred2.graph <- ggplot(data=test) +
geom_rect(aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470)) +
scale_fill_manual(values = cols)
pred2.graph
pred2 <- factor(pred, levels = c("none", "few", "few more", "little", "some", "most", "all"))
cols = c("none"="#ceb467", "few" = "#ace5b2", "few more" = "#7aef87", "little" = "#57d165", "some" = "#21a31d", "most" = "#197f24", "all" = "#115118")
pred2.graph <- ggplot(data=test) +
geom_rect(aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470)) +
scale_fill_manual(values = cols)
pred2.graph
species.map$scam <- factor(species.map$scam, levels = c("none", "few", "few more", "little", "some", "most", "all"))
cols = c("none"="#ceb467", "few" = "#ace5b2", "few more" = "#7aef87", "little" = "#57d165", "some" = "#21a31d", "most" = "#197f24", "all" = "#115118")
serc.plots <- ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Population Abundance of SERC Plots"), x="Easting", y="Northing", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470)) +
scale_fill_manual(values = cols)
serc.plots
grid.arrange(serc.plots, pred2.graph, ncol=2)
ggsave("plots/plot_comparison.png", arrangeGrob(serc.plots, pred2.graph, ncol=2), width = 11)
serc.plots <- ggplot() +
geom_rect(data=species.map, aes(xmin=easting, xmax=easting + 20, ymin=northing, ymax=northing + 20, fill = as.factor(unlist(species.map[1]))), color=NA) +
labs(title=paste("Scam Population Abundance in SERC Plots"), x="Easting", y="Northing", fill = "Cover") +
geom_rect(data=training, aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30), color="black", fill = NA) +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470)) +
scale_fill_manual(values = cols)
serc.plots
ggsave("plots/prediction_validation.png")
pred2.graph <- ggplot(data=test) +
geom_rect(aes(xmin=easting, xmax=easting + 30, ymin=northing, ymax=northing + 30, fill = as.factor(unlist(pred2))), color="black") +
labs(title=paste(curr.species, "Population Abundance Predictions"), x="Easting", y="Northing", fill="Cover") +
scale_x_continuous(limits = c(365430, 366280)) +
scale_y_continuous(limits = c(4303800, 4304470)) +
scale_fill_manual(values = cols)
pred2.graph
ggsave("plots/prediction_visual2.png")
# need to run training.R concurrently
grid.arrange(serc.plots, pred2.graph, ncol=2)
ggsave("plots/plot_comparison.png", arrangeGrob(serc.plots, pred2.graph, ncol=2), width = 11)
ggplot(data = full.data, aes(x=x, y=y)) + geom_point(aes(color = ndvi))
ggsave("plots/score_test.png")
ggsave("plots/score_test.png")
ggsave("plots/ndvi_test.png")
int("1")
as.integer("1")
load(lake.inputs)
lake.inputs = paste("InputsCV", lake.name, ".rda",sep="")
lake.name <- chopper
lake.name <- "chopper"
lake.inputs = paste("InputsCV", lake.name, ".rda",sep="")
load(lake.inputs)
setwd("~/Documents/Junior_Year/DISC_REU/DISC_bayesian_model/")
lake.inputs = paste("InputsCV", lake.name, ".rda",sep="")
load(lake.inputs)
run.cv = function(x){
attach(x)
out = try(adaptMetropGibbs(ltd=ltd,starting=starting,tuning=tune,batch=2400,report=400,
inputs=inputs,priors=priors))
n.sim = 2400*25
n.burn = 50000
n.step = 5
save.idx = seq(n.burn+1,n.sim,n.step)
if (class(out)=="try-error"){
samps = array(NA,dim=c(length(save.idx),length(starting)))
} else {
samps = out$p.theta.samples[save.idx,]
}
return(samps)
detach(x)
}
####### END OF Malcolm's CODE ######
results = vector("list", ((stop - start) + 1))
for (i in start+1 : stop){
results[[i]] = run.cv(mod.inputs[i])
}
return(results)
start <- 1
stop <- 3
run.cv = function(x){
attach(x)
out = try(adaptMetropGibbs(ltd=ltd,starting=starting,tuning=tune,batch=2400,report=400,
inputs=inputs,priors=priors))
n.sim = 2400*25
n.burn = 50000
n.step = 5
save.idx = seq(n.burn+1,n.sim,n.step)
if (class(out)=="try-error"){
samps = array(NA,dim=c(length(save.idx),length(starting)))
} else {
samps = out$p.theta.samples[save.idx,]
}
return(samps)
detach(x)
}
####### END OF Malcolm's CODE ######
results = vector("list", ((stop - start) + 1))
for (i in start+1 : stop){
results[[i]] = run.cv(mod.inputs[i])
}
return(results)
library(spBayes)
library(mgcv)
run.cv = function(x){
attach(x)
out = try(adaptMetropGibbs(ltd=ltd,starting=starting,tuning=tune,batch=2400,report=400,
inputs=inputs,priors=priors))
n.sim = 2400*25
n.burn = 50000
n.step = 5
save.idx = seq(n.burn+1,n.sim,n.step)
if (class(out)=="try-error"){
samps = array(NA,dim=c(length(save.idx),length(starting)))
} else {
samps = out$p.theta.samples[save.idx,]
}
return(samps)
detach(x)
}
####### END OF Malcolm's CODE ######
results = vector("list", ((stop - start) + 1))
for (i in start+1 : stop){
results[[i]] = run.cv(mod.inputs[i])
}
return(results)
results = vector("list", ((stop - start) + 1))
for (i in start+1 : stop){
results[[i]] = run.cv(mod.inputs[[i]])
}
start+1
for (i in start : stop){
results[[i]] = run.cv(mod.inputs[[i]])
}
start <- 0
for (i in start+1 : stop){
results[[i]] = run.cv(mod.inputs[[i]])
}
results = vector("list", ((stop - start) + 1))
for (i in start+1 : stop){
results[[i]] = run.cv(mod.inputs[[i]])
}
View(training)
min(training$overlap)
rm(list=ls())
# import libraries
library(spBayes)
library(ggplot2)
library(coda)
library(gridExtra)
library(mgcv)
library(zoo)
library(Matrix)
library(stringr)
library(snow)
library(rlecuyer)
library(snowfall)
library(readxl)
#STEP 1: Setup and Data Manipulation
setwd("~/Documents/Junior_Year/DISC_REU/DISC_bayesian_model/")
# specify sheet number
curr.lake <- read_excel("bigwoods.xls", sheet = 5, col_names = TRUE)
char.dat <- data.frame(age=rep(NA,length(curr.lake$Date)))
# define variables from lake
char.dat$age <- with(curr.lake, round(Date) - 2018) # calculate YBP
char.dat$sed.rate <- with(curr.lake, Flux / Count) # NOTE: do not have age of sediment core so we left it fixed
char.dat$influx <- curr.lake$`Char Flux`
char.dat$count <- curr.lake$Count
qplot(age,count,data=char.dat)  + geom_smooth(method="loess", span=0.08) + theme_bw()
# specify sheet number
curr.lake <- read_excel("bigwoods.xls", sheet = 7, col_names = TRUE)
char.dat <- data.frame(age=rep(NA,length(curr.lake$Date)))
# define variables from lake
char.dat$age <- with(curr.lake, round(Date) - 2018) # calculate YBP
char.dat$sed.rate <- with(curr.lake, Flux / Count) # NOTE: do not have age of sediment core so we left it fixed
char.dat$influx <- curr.lake$`Char Flux`
char.dat$count <- curr.lake$Count
qplot(age,count,data=char.dat)  + geom_smooth(method="loess", span=0.08) + theme_bw()
# specify sheet number
curr.lake <- read_excel("bigwoods.xls", sheet = 5, col_names = TRUE)
char.dat <- data.frame(age=rep(NA,length(curr.lake$Date)))
# define variables from lake
char.dat$age <- with(curr.lake, round(Date) - 2018) # calculate YBP
char.dat$sed.rate <- with(curr.lake, Flux / Count) # NOTE: do not have age of sediment core so we left it fixed
char.dat$influx <- curr.lake$`Char Flux`
char.dat$count <- curr.lake$Count
qplot(age,count,data=char.dat)  + geom_smooth(method="loess", span=0.08) + theme_bw()
#STEP 2: Approximate background and foreground intensity (using cubic base spline) and set offset
char.dat$age.c = with(char.dat,age-min(age))
char.dat$age.s = with(char.dat,age.c/max(age.c))
n = nrow(char.dat)
n.knots = 31 # this variable will change for each model based on time interval
# note: smoothCon is part of the mgcv package for constructing the smooth terms in a GAM model
CRbasis = smoothCon(s(age.s,k=n.knots,bs="cr"),data=char.dat,knots=NULL,absorb.cons=TRUE,
scale.penalty=TRUE)
Sb = CRbasis[[1]]$S[[1]]
X = CRbasis[[1]]$X
knots = as.numeric(CRbasis[[1]]$xp)
S.scale = CRbasis[[1]]$S.scale
TT <- 1/lake.dat$sed.rate # sets offset (inverse of sedimentation rate from STEP 1)
char.dat$t.hold = rollapply(char.dat$count,50,function(x)quantile(x,0.9),fill="extend")
# Separate charcoal counts
y.back = y.fore = char.dat$count
y.back[char.dat$count>char.dat$t.hold] = NA # removes count value for background if the value is distributed in the 90% quantile
y.fore[char.dat$count<char.dat$t.hold] = NA # removes count value for background if the value is NOT distributed in the 90% quantile
#approximate the missing values (background = na.approx, foreground = randomly pulls value from binomial distribution)
y.back = round(na.approx(y.back))
y.fore[is.na(y.fore)] = rbinom(sum(is.na(y.fore)),1,0.1) # ask about why we add a bunch of zeros instead of NAs
# Plot background and foreground counts
plot(char.dat$age,y.back)
# Plot background and foreground counts
plot(char.dat$age,y.back)
# Plot background and foreground counts
plot(char.dat$age,y.back)
plot(char.dat$age,y.fore)
# Add separate counts to data.frame
char.dat$count.b = y.back
char.dat$count.f = y.fore
# Fit GAM Model to background (gam is part of mgcv package)
m1 = gam(y.back~X,family=poisson,offset=log(TT),paraPen=list(X=list(Sb,sp=0.001*S.scale)))
# plots the points and our fitted GAM model for background
plot(char.dat$age,y.back)
lines(char.dat$age,fitted(m1),col="blue",lwd=2)
# Fit GAM Model to foreground
m2 = gam(y.fore~X,family=poisson,offset=log(TT),paraPen=list(X=list(Sb,sp=(1e-7)*S.scale)))
# plots the points and our fitted GAM model for foreground
plot(char.dat$age,y.fore)
lines(char.dat$age,fitted(m2),col="blue",lwd=2)
rm(list=ls())