library(mapproj)
library(ggplot2)
install.packages("map")
library(map)
install.packages("MAP")
library(MAP)
install.packages("cowplot")
install.packages("cowplot")
library(cowplot)
install.packages("mapproj")
install.packages("mapproj")
library(mapproj)
library(ggplot2)
# retrieve map data for USA
us_map <- map_data("usa")
# create plot for map
p_map <- ggplot() +
geom_polygon(data = us_map, aes(x = long, y = lat, group = group)) +
coord_map()
ggplot() +
geom_polygon(data = us_map, aes(x = long, y = lat, group = group)) +
coord_map()
gplot(data, aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "red")
ggplot(data, aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "red")
p_data <- ggplot(data, aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "red")
combined_plot <- p_map +
geom_tile(data = data, aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
scale_fill_gradient(low = "white", high = "red")
# display combined plot
print(combined_plot)
p_map +
geom_tile(data = data, aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
scale_fill_gradient(low = "white", high = "red")
p_map <- ggplot() +
geom_polygon(data = us_map, aes(x = long, y = lat, group = group)) +
coord_map()
# create plot for data
p_data <- ggplot(data, aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "red")
# combine plots
combined_plot <- p_map +
geom_tile(data = data, aes(x = LON, y = LAT, fill = p_QNITROgl)) +
scale_fill_gradient(low = "white", high = "red")
# display combined plot
print(combined_plot)
p_map <- ggplot() +
geom_polygon(data = us_map, aes(x = long, y = lat, group = group)) +
coord_map()
# create plot for data
p_data <- ggplot(data, aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "blue")
# combine plots
combined_plot <- p_map +
geom_tile(data = data, aes(x = LON, y = LAT, fill = p_QNITROgl)) +
scale_fill_gradient(low = "white", high = "blue")
# display combined plot
print(combined_plot)
p_data <- ggplot(data = subset(data, LON >= -60), aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "red") +
coord_cartesian(xlim = c(-60, NA))
ggplot(data = subset(data, LON >= -60), aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
+     geom_tile() +
+     scale_fill_gradient(low = "white", high = "red") +
+     coord_cartesian(xlim = c(-60, NA))
# read data from Excel file
data <- read_excel("C:\\Mini-Model-TEST\\results2.xlsx")
# retrieve map data for USA
us_map <- map_data("usa")
# create map
ggplot() +
geom_polygon(data = us_map, aes(x = long, y = lat, group = group))
install.packages("MAP")
install.packages("MAP")
library(MAP)
install.packages("cowplot")
install.packages("cowplot")
library(cowplot)
install.packages("mapproj")
install.packages("mapproj")
library(mapproj)
library(ggplot2)
# retrieve map data for USA
us_map <- map_data("usa")
# create plot for map
p_map <- ggplot() +
geom_polygon(data = us_map, aes(x = long, y = lat, group = group)) +
coord_map()
# create plot for data
p_data <- ggplot(data, aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "red")
# combine plots
combined_plot <- p_map +
geom_tile(data = data, aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
scale_fill_gradient(low = "white", high = "red")
# display combined plot
print(combined_plot)
p_map <- ggplot() +
geom_polygon(data = us_map, aes(x = long, y = lat, group = group)) +
coord_map()
# create plot for data
p_data <- ggplot(data, aes(x = LON, y = LAT, fill = p_QNLEACHgl)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "blue")
# combine plots
combined_plot <- p_map +
geom_tile(data = data, aes(x = LON, y = LAT, fill = p_QNITROgl)) +
scale_fill_gradient(low = "white", high = "blue")
# display combined plot
print(combined_plot)
library(MASS)
set.seed(123)##sets random
beta_true <- c(2, 0.5)
phi_true <- 2.5
num_simulations <- 250 ##number of lines generates
###################
#####Round One#####
###################
n <- 50
x1 <- rnorm(n, mean=2,sd = sqrt(1/phi_true))##sets x1 norm around 2
X <- cbind(1, x1)## creates X matrix
mu_true <- exp(beta_true[1] + beta_true[2] * x1)
y<-rnegbin(n, mu = mu_true, theta = 2.5)
start_params <- c(2, 0.5, log(2.5))
negbin_likelihood <- function(params, y, X) {
beta <-params[1:2]
r <- exp(params[3])
mu <- exp(X %*% beta)
size_factor <- r/(r + mu)
loglik <- sum(lgamma(y + r) - lgamma(y + 1) - lgamma(r) + r*log(size_factor) + y*log(1-size_factor))
return(-loglik)
}
results1 <- data.frame(matrix(NA, nrow = num_simulations, ncol = length(start_params)))
for (i in 1:num_simulations) {
x1 <- rnorm(n,mean=2,sd = sqrt(1/phi_true))
mu_true <- exp(beta_true[1] + beta_true[2]*x1)
y<-rnegbin(n, mu = mu_true, theta = 2.5)
data_sim <- data.frame(x1 = x1, y = y)
start_params <-  c(2, 0.5, log(2.5))
gfit <- optim(start_params, negbin_likelihood, y = y, X = cbind(1, x1), method = "BFGS")
results1[i, 1] <- gfit$par[1]
results1[i, 2] <- gfit$par[2]
results1[i, 3] <- exp(gfit$par[3])
}
# Calculate the mean and standard deviation of the estimated coefficients across all simulations
coef_mean <- colMeans(results1)
coef_sd <- apply(results1, 2, sd)
# Print the results
cat("True coefficients:", paste(beta_true, collapse = ", "))
cat("Estimated coefficients (mean):", paste(round(coef_mean, 3), collapse = ", "))
cat("Standard deviation of estimated coefficients:", paste(round(coef_sd, 3), collapse = ", "))
###################
#####Round Two#####
###################
n <- 1000
x1 <- rnorm(n, mean=2,sd = sqrt(1/phi_true))##sets x1 norm around 2
X <- cbind(1, x1)## creates X matrix
mu_true <- exp(beta_true[1] + beta_true[2] * x1)
y<-rnegbin(n, mu = mu_true, theta = 2.5)
start_params <- c(2, 0.5, log(2.5))
results2 <- data.frame(matrix(NA, nrow = num_simulations, ncol = length(start_params)))
for (i in 1:num_simulations) {
x1 <- rnorm(n,mean=2,sd = sqrt(1/phi_true))
mu_true <- exp(beta_true[1] + beta_true[2]*x1)
y<-rnegbin(n, mu = mu_true, theta = 2.5)
data_sim <- data.frame(x1 = x1, y = y)
start_params <-  c(2, 0.5, log(2.5))
gfit <- optim(start_params, negbin_likelihood, y = y, X = cbind(1, x1), method = "BFGS")
results2[i, 1] <- gfit$par[1]
results2[i, 2] <- gfit$par[2]
results2[i, 3] <- exp(gfit$par[3])
}
# Calculate the mean and standard deviation of the estimated coefficients across all simulations
coef_mean <- colMeans(results2)
coef_sd <- apply(results2, 2, sd)
# Print the results
cat("True coefficients:", paste(beta_true, collapse = ", "))
cat("Estimated coefficients (mean):", paste(round(coef_mean, 3), collapse = ", "))
cat("Standard deviation of estimated coefficients:", paste(round(coef_sd, 3), collapse = ", "))
#####################
#####Round Three#####
#####################
n <- 2000
x1 <- rnorm(n, mean=2,sd = sqrt(1/phi_true))##sets x1 norm around 2
X <- cbind(1, x1)## creates X matrix
mu_true <- exp(beta_true[1] + beta_true[2] * x1)
y<-rnegbin(n, mu = mu_true, theta = 2.5)
start_params <- c(2, 0.5, log(2.5))
results3 <- data.frame(matrix(NA, nrow = num_simulations, ncol = length(start_params)))
for (i in 1:num_simulations) {
x1 <- rnorm(n,mean=2,sd = sqrt(1/phi_true))
mu_true <- exp(beta_true[1] + beta_true[2]*x1)
y<-rnegbin(n, mu = mu_true, theta = 2.5)
data_sim <- data.frame(x1 = x1, y = y)
gfit <- optim(start_params, negbin_likelihood, y = y, X = cbind(1, x1), method = "BFGS")
results3[i, 1] <- gfit$par[1]
results3[i, 2] <- gfit$par[2]
results3[i, 3] <- exp(gfit$par[3])
}
# Calculate the mean and standard deviation of the estimated coefficients across all simulations
coef_mean <- colMeans(results3)
coef_sd <- apply(results3, 2, sd)
# Print the results
cat("True coefficients:", paste(beta_true, collapse = ", "))
cat("Estimated coefficients (mean):", paste(round(coef_mean, 3), collapse = ", "))
cat("Standard deviation of estimated coefficients:", paste(round(coef_sd, 3), collapse = ", "))
###################
#####Graphs#####
###################
library(ggplot2)
library(dplyr)
# Convert the results matrix to a data frame
results1_df <- data.frame(results1) %>% rename(B1 = X1, B2 = X2, size = X3)%>% select(-size)
results2_df <- data.frame(results2) %>% rename(B1 = X1, B2 = X2, size = X3)%>% select(-size)
results3_df <- data.frame(results3) %>% rename(B1 = X1, B2 = X2, size = X3)%>% select(-size)
# Reshape the data from wide to long format
results1_long <- tidyr::pivot_longer(results1_df, cols = 1:2, names_to = "param", values_to = "estimate")
results2_long <- tidyr::pivot_longer(results2_df, cols = 1:2, names_to = "param", values_to = "estimate")
results3_long <- tidyr::pivot_longer(results3_df, cols = 1:2, names_to = "param", values_to = "estimate")
# Create four individual density plots of the estimated coefficients
ggplot(results1_long, aes(x = estimate, color = "Simulation 1")) +
geom_density(alpha = 0.6) +
labs(x = "Estimated Coefficient", y = "Density", color = "Simulation") +
facet_wrap(~param, ncol = 2, scales = "free") +
geom_density(data = results2_long, aes(x = estimate, color = "Simulation 2"), alpha = 0.6) +
geom_density(data = results3_long, aes(x = estimate, color = "Simulation 3"), alpha = 0.6)
# Filter results to only include B1 parameter
results1_long_B1 <- results1_long %>% filter(param == "B1")
results2_long_B1 <- results2_long %>% filter(param == "B1")
results3_long_B1 <- results3_long %>% filter(param == "B1")
# Create density plot for B1 across all simulations
ggplot(results1_long_B1, aes(x = estimate, color = "Simulation 1")) +
geom_density(alpha = 0.6) +
labs(x = "Estimated Coefficient (B1)", y = "Density", color = "Simulation") +
geom_density(data = results2_long_B1, aes(x = estimate, color = "Simulation 2"), alpha = 0.6) +
geom_density(data = results3_long_B1, aes(x = estimate, color = "Simulation 3"), alpha = 0.6) +
scale_color_manual(values = c("Simulation 1" = "red", "Simulation 2" = "blue", "Simulation 3" = "green"))+
geom_vline(xintercept = beta_true[1], linetype = "dashed", color = "black")
# Filter results to only include B2 parameter
results1_long_B2 <- results1_long %>% filter(param == "B2")
results2_long_B2 <- results2_long %>% filter(param == "B2")
results3_long_B2 <- results3_long %>% filter(param == "B2")
# Create density plot for B2 across all simulations
ggplot(results1_long_B2, aes(x = estimate, color = "Simulation 1")) +
geom_density(alpha = 0.6) +
labs(x = "Estimated Coefficient (B2)", y = "Density", color = "Simulation") +
geom_density(data = results2_long_B2, aes(x = estimate, color = "Simulation 2"), alpha = 0.6) +
geom_density(data = results3_long_B2, aes(x = estimate, color = "Simulation 3"), alpha = 0.6) +
scale_color_manual(values = c("Simulation 1" = "red", "Simulation 2" = "blue", "Simulation 3" = "green")) +
geom_vline(xintercept = beta_true[2], linetype = "dashed", color = "black")
# ---- 0. Libraries, \in files, setwd " ----
rm(list = ls())
setwd("C:/Mini-Model-Redo/R")
library(rgdal) #redOGR,
library(raster) #intersect
library(sp) #intersect
library(readxl) # read_excel
library(writexl) # write_excel
library(rgdal) #redOGR,
library(raster) #intersect
datg.list   <- list()
mar_def     <- c(5.1, 4.1, 4.1, 2.1) # default margins
state_us      <- readOGR("in\\cb_2018_us_state_500k\\cb_2018_us_state_500k.shp")
newcrs        <- CRS("+proj=longlat +datum=WGS84")
state_us      <- spTransform(state_us, newcrs)
e <- extent(-125, -59, 25, 50)
state_cous <- crop(state_us, e)
SNO_GID       <- read_excel("in\\SNO_GRID.xlsx")
View(SNO_GID)
# ---- 0. Libraries, \in files, setwd " ----
rm(list = ls())
setwd("C:/Mini-Model-Redo/R")
library(rgdal) #redOGR,
library(raster) #intersect
library(sp) #intersect
library(readxl) # read_excel
library(writexl) # write_excel
library(rgdal) #redOGR,
library(raster) #intersect
datg.list   <- list()
mar_def     <- c(5.1, 4.1, 4.1, 2.1) # default margins
state_us      <- readOGR("in\\cb_2018_us_state_500k\\cb_2018_us_state_500k.shp")
newcrs        <- CRS("+proj=longlat +datum=WGS84")
# ---- 0. Libraries, \in files, setwd " ----
rm(list = ls())
setwd("C:/Mini-Model-Redo/R")
library(rgdal) #redOGR,
library(raster) #intersect
library(sp) #intersect
library(readxl) # read_excel
library(writexl) # write_excel
library(rgdal) #redOGR,
library(raster) #intersect
datg.list   <- list()
mar_def     <- c(5.1, 4.1, 4.1, 2.1) # default margins
state_us      <- readOGR("in\\cb_2018_us_state_500k\\cb_2018_us_state_500k.shp")
newcrs        <- CRS("+proj=longlat +datum=WGS84")
state_us      <- spTransform(state_us, newcrs)
e <- extent(-125, -59, 25, 50)
state_cous <- crop(state_us, e)
SNO_GID       <- read_excel("in\\SNO_GRID.xlsx")
# ---- 1.  Define experiment name  [[1]]  "simpleg_maxi"----
cmf <- as.character("simpleg_maxi") # name of cmf file
out <- as.character("out\\maxi") # name of out folder
# --- Run GEMPACK Results2Excel.R file ---
source("ResultsGCZ2Excel.R")
# ---- 0. Libraries, \in files, setwd " ----
rm(list = ls())
setwd("C:/Replication Chapter/SIMPLE_G1")
library(rgdal) #redOGR,
library(raster) #intersect
library(sp) #intersect
library(readxl) # read_excel
library(writexl) # write_excel
library(rgdal) #redOGR,
library(raster) #intersect
datg.list   <- list()
mar_def     <- c(5.1, 4.1, 4.1, 2.1) # default margins
state_us      <- readOGR("in\\cb_2018_us_state_500k\\cb_2018_us_state_500k.shp")
newcrs        <- CRS("+proj=longlat +datum=WGS84")
state_us      <- spTransform(state_us, newcrs)
# ---- 0. Libraries, \in files, setwd " ----
rm(list = ls())
setwd("C:/Replication Chapter/SIMPLE_G1/R")
library(rgdal) #redOGR,
library(raster) #intersect
library(sp) #intersect
library(readxl) # read_excel
library(writexl) # write_excel
library(rgdal) #redOGR,
library(raster) #intersect
datg.list   <- list()
mar_def     <- c(5.1, 4.1, 4.1, 2.1) # default margins
state_us      <- readOGR("in\\cb_2018_us_state_500k\\cb_2018_us_state_500k.shp")
newcrs        <- CRS("+proj=longlat +datum=WGS84")
state_us      <- spTransform(state_us, newcrs)
e <- extent(-125, -59, 25, 50)
state_cous <- crop(state_us, e)
SNO_GID       <- read_excel("in\\SNO_GRID.xlsx")
# ---- 1.  Define experiment name  [[1]]  "simpleg_maxi"----
cmf <- as.character("simpleg_maxi") # name of cmf file
out <- as.character("out\\maxi") # name of out folder
# --- Run GEMPACK Results2Excel.R file ---
source("ResultsGCZ2Excel.R")
# Load required libraries if not already loaded
library(rgdal)
library(raster)
library(sp)
library(readxl)
library(writexl)
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the data from the CSV file
qcrop_data <- read.csv("qcrop.cvs")
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the data from the CSV file
qcrop_data <- read.csv("qcrop.cvs")
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the data from the CSV file
qcrop_data <- read.csv("qcrop.cvs")
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the data from the CSV file
qcrop_data <- read.csv("qcrop.cvs")
# Create a SpatialPointsDataFrame from the data
coordinates(qcrop_data) <- c("long", "lat")
# Load required libraries if not already loaded
library(rgdal)
library(raster)
library(sp)
library(readxl)
library(writexl)
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the data from the CSV file
qcrop_data <- read.csv("qcrop.cvs")
# Read the data from the CSV file
qcrop_data <- read.csv("C:/Replication Chapter/SIMPLE_G1/R/qcrop.csv")
# Create a SpatialPointsDataFrame from the data
coordinates(qcrop_data) <- c("long", "lat")
str(qcrop_data)
library(writexl)
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the data from the CSV file
qcrop_data <- read.csv("C:/Replication Chapter/SIMPLE_G1/R/qcrop.csv")
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the data from the CSV file
qcrop_data <- read.csv("C:/Replication Chapter/SIMPLE_G1/R/qcrop.csv")
str(qcrop_data)
# Create a SpatialPointsDataFrame from the data
coordinates(qcrop_data) <- c("long", "lat")
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the data from the CSV file
qcrop_data <- read.csv("C:/Replication Chapter/SIMPLE_G1/R/qcrop.csv")
str(qcrop_data)
# Create a SpatialPointsDataFrame from the data
coordinates(qcrop_data) <- c("long", "lat")
# Define the CRS for the SpatialPointsDataFrame
proj4string(qcrop_data) <- CRS("+proj=longlat +datum=WGS84")
# Create the plot
png(filename = "plot\\QCROP_QLAND.png", units = "in", width = 20, height = 10, res = 300)
par(mfrow = c(1, 2), mar = mar_def)
plot(qcrop_data$QCROP, main = "p_QCROP_i", col = qcrop2i.col, breaks = qcrop2i.brk,
box = FALSE, axes = FALSE, legend = FALSE)
# Load necessary libraries if not loaded
library(sp)
library(raster)
# Load necessary libraries if not loaded
library(sp)
library(raster)
# Read the CSV file
qcrop_data <- read.csv("qcrop.csv")
# Merge qnew_variable data with GID data (assuming 'SNO' is the common identifier)
merged_data <- merge(qcrop_data, SNO_GID, by = "SNO")
# Convert to spatial points
coordinates(merged_data) <- c("Longitude", "Latitude")
# Load necessary libraries if not loaded
library(sp)
library(raster)
library(readxl)
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the CSV file
qcrop_data <- read.csv("qcrop.csv")
# Convert qcrop_data to a spatial points data frame
coordinates(qcrop_data) <- c("long", "lat")
# Create a raster stack
qnew_raster <- stack(qcrop_data)
# Check unique values in Longitude and Latitude columns
unique(qcrop_data$Longitude_column_name)
unique(qcrop_data$Latitude_column_name)
# Check column names in qcrop_data
names(qcrop_data)
View(qcrop_data)
View(qcrop_data)
View(qcrop_data)
# Load necessary libraries if not loaded
library(sp)
library(raster)
library(readxl)
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the CSV file
qcrop_data <- read.csv(""C:/Replication Chapter/SIMPLE_G1/R/qcrop.csv"")
# Read the CSV file
qcrop_data <- read.csv("C:/Replication Chapter/SIMPLE_G1/R/qcrop.csv")
# Convert qcrop_data to a spatial points data frame
coordinates(qcrop_data) <- c("long", "lat")
# Create a raster stack
qnew_raster <- stack(qcrop_data)
unique(qcrop_data$qcrop)
library(raster)
# Convert 'qcrop' column to numeric
qcrop_data$qcrop <- as.numeric(as.character(qcrop_data$qcrop))
# Create a raster stack using the 'qcrop' column
qcrop_raster <- rasterFromXYZ(qcrop_data[, c("Longitude", "Latitude", "qcrop")])
# Create a raster stack using the 'qcrop' column
qcrop_raster <- rasterFromXYZ(qcrop_data[, c("long", "lat", "qcrop")])
View(qcrop_data)
# Load necessary libraries if not loaded
library(sp)
library(raster)
library(readxl)
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the CSV file
qcrop_data <- read.csv("C:/Replication Chapter/SIMPLE_G1/R/qcrop.csv")
# Convert qcrop_data to a spatial points data frame
coordinates(qcrop_data) <- c("long", "lat")
# Load necessary libraries if not loaded
library(sp)
library(raster)
library(readxl)
# Set the working directory
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the CSV file
qcrop_data <- read.csv("qcrop.csv")
# Convert qcrop_data to a spatial points data frame
coordinates(qcrop_data) <- c("long", "lat")
# Load necessary libraries if not loaded
library(sp)
library(raster)
library(readxl)
# Set the working directory
setwd("C:/Replication Chapter/SIMPLE_G1/R")
# Read the CSV file
qcrop_data <- read.csv("qcrop.csv")
# Convert 'qcrop' column to numeric
qcrop_data$qcrop <- as.numeric(as.character(qcrop_data$qcrop))
# Create a SpatialPointsDataFrame
coordinates(qcrop_data) <- c("long", "lat")  # Set the coordinates
library(sp)
library(raster)
library(readxl)
setwd("C:/Replication Chapter/SIMPLE_G1/R")
qcrop_data <- read.csv("qcrop.csv")
# Convert 'qcrop' column to numeric
qcrop_data$qcrop <- as.numeric(as.character(qcrop_data$qcrop))
# Create SpatialPointsDataFrame
qcrop_sp <- SpatialPointsDataFrame(
coords = qcrop_data[, c("long", "lat")],
data = qcrop_data,
proj4string = CRS("+proj=longlat +datum=WGS84")  # Define the projection
)
View(qcrop_sp)
# Create a raster stack using the 'qcrop' column
qcrop_raster <- rasterFromXYZ(qcrop_sp)
# Save the raster stack to a file
writeRaster(qcrop_raster, filename = "qcrop_raster.tif", format = "GTiff")
