This repository has been archived by the owner on Jun 13, 2023. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 18
EBS Jellyfish
Jim Thorson edited this page Oct 25, 2018
·
3 revisions
- Note: Ellen Yasumiishi added this, James Thorson is mildly editing to improve readability.
setwd("C:/Users/Alien/Documents")
temp <- read.csv("EBSJellyfish.csv")
Sci=temp[,1] #Species name
year = temp[,2]
vessel = temp[,3] #Set to 1
AreaSwept = temp[,4]
Lat = temp[,5] #Latitude start or equilibrium in decimal degrees
Lon = temp[,6]
Wt=temp[,7] #Kilograms
Version = "VAST_v5_1_0"
Method = c("Grid", "Mesh")[2]
grid_size_km = 30
n_x = c(50, 100, 200, 500, 1000, 2000)[1]
Kmeans_Config = list( "randomseed"=1, "nstart"=100, "iter.max"=1e3 )
FieldConfig = c("Omega1"=5, "Epsilon1"=5, "Omega2"=5, "Epsilon2"=5)
RhoConfig = c("Beta1"=0, "Beta2"=0, "Epsilon1"=0, "Epsilon2"=0)
OverdispersionConfig = c("Vessel"=0, "VesselYear"=0)
ObsModel = c(1,1)
Options = c(SD_site_density = 0, SD_site_logdensity = 0,
Calculate_Range = 1, Calculate_evenness = 0,
Calculate_effective_area = 1, Calculate_Cov_SE = 0,
Calculate_Synchrony = 1, Calculate_Coherence = 1)
strata.limits <- data.frame(STRATA = "EBS")
Data_Set=temp
Species_set=c("Aequorea", "Aurelia", "Chrysaora", "Cyanea", "Staurophora")
library(pander)
pander::pandoc.table(summary(Sci))
DateFile = paste0(getwd(),'/VAST_Jellyfish2017/')
dir.create(DateFile)
Record = list(Version = Version, Method = Method,
grid_size_km = grid_size_km, n_x = n_x,
FieldConfig = FieldConfig, RhoConfig = RhoConfig,
OverdispersionConfig = OverdispersionConfig,
ObsModel = ObsModel, Kmeans_Config = Kmeans_Config,
Region = Region, Species_set = Species_set, strata.limits = strata.limits)
save(Record, file = file.path(DateFile, "Record.RData"))
capture.output(Record, file = paste0(DateFile, "Record.txt"))
Data_Geostat = data.frame(spp = temp[, "Sci"], Year = temp[, "year"], Catch_KG = temp[, "Wt"], AreaSwept_km2= temp[, "AreaSwept"], Vessel= temp[, "vessel"], Lat = temp[, "Lat"], Lon = temp[, "Lon"])
pander::pandoc.table( head(Data_Geostat), digits=3 )
if (Region == "Other") { Extrapolation_List = SpatialDeltaGLMM::Prepare_Extrapolation_Data_Fn(Region = Region,
strata.limits = strata.limits, observations_LL = Data_Geostat[, c("Lat", "Lon")], maximum_distance_from_sample = 37.5) }
Spatial_List = SpatialDeltaGLMM::Spatial_Information_Fn( grid_size_km=grid_size_km,
n_x=n_x, Method=Method, Lon=Data_Geostat[,"Lon"],
Lat=Data_Geostat[,"Lat"], Extrapolation_List=Extrapolation_List,
randomseed=Kmeans_Config[["randomseed"]], nstart=Kmeans_Config[["nstart"]],
iter.max=Kmeans_Config[["iter.max"]], DirPath=DateFile,
Save_Results=FALSE)
Data_Geostat = cbind( Data_Geostat, knot_i=Spatial_List$knot_i)
TmbData = Data_Fn(Version=Version, FieldConfig=FieldConfig,RhoConfig=RhoConfig, ObsModel=ObsModel,
b_i=Data_Geostat[,'Catch_KG'], a_i=Data_Geostat[,'AreaSwept_km2'],c_i=as.numeric(Data_Geostat[,'spp'] )-1,
v_i=as.numeric(Data_Geostat[,'Vessel'])-1, s_i=Data_Geostat[,'knot_i']-1,
t_i=Data_Geostat[,'Year'], a_xl=Spatial_List$a_xl,
MeshList=Spatial_List$MeshList, GridList=Spatial_List$GridList,
Method=Spatial_List$Method, Options=Options )
TmbList = Build_TMB_Fn(TmbData=TmbData, RunDir=DateFile,
Version=Version, RhoConfig=RhoConfig, loc_x=Spatial_List$loc_x,Method=Method)
Obj = TmbList[["Obj"]]
SpatialDeltaGLMM::Plot_data_and_knots(Extrapolation_List = Extrapolation_List,
Spatial_List=Spatial_List,Data_Geostat=Data_Geostat,
PlotDir=DateFile)
Opt = TMBhelper::Optimize( obj=Obj, lower=TmbList[["Lower"]],
upper=TmbList[["Upper"]], getsd=TRUE, savedir=DateFile,
bias.correct=TRUE, newtonsteps=1 )
Report = Obj$report()
Save = list("Opt"=Opt, "Report"=Report, "ParHat"=Obj$env$parList(Opt$par), "TmbData"=TmbData)
save(Save, file=paste0(DateFile,"Save.RData"))
library(pander)
pander::pandoc.table(Opt$diagnostics[,c('Param','Lower','MLE','Upper','final_gradient')])
Enc_prob=SpatialDeltaGLMM::Check_encounter_prob(Report=Report,
Data_Geostat=Data_Geostat,DirName=DateFile)
Q=SpatialDeltaGLMM::QQ_Fn(TmbData=TmbData, Report=Report,
FileName_PP=paste0(DateFile,"Posterior-Predictive.jpg"),
FileName_Phist=paste0(DateFile,"Posterior_Predictive-Histogram.jpg"),
FileName_QQ=paste0(DateFile,"Q-Q_plot.jpg"),
FileName_Qhist=paste0(DateFile,"Q-Q_hist.jpg")) #SpatialDeltaGLMM::
MapDetails_List = SpatialDeltaGLMM::MapDetails_Fn( "Region"=Region, "NN_Extrap"=Spatial_List$PolygonList$NN_Extrap, "Extrapolation_List"=Extrapolation_List )
Year_Set = seq(min(Data_Geostat[,'Year']),max(Data_Geostat[,'Year']))
Years2Include = which( Year_Set %in% sort(unique(Data_Geostat[,'Year'])))
SpatialDeltaGLMM::plot_residuals(Lat_i=Data_Geostat[,
"Lat"], Lon_i=Data_Geostat[,"Lon"],TmbData=TmbData,
Report=Report, Q=Q, savedir = DateFile, MappingDetails=MapDetails_List[["MappingDetails"]],
PlotDF=MapDetails_List[["PlotDF"]],MapSizeRatio=MapDetails_List[["MapSizeRatio"]],
Xlim=MapDetails_List[["Xlim"]],Ylim=MapDetails_List[["Ylim"]],
FileName=DateFile,Year_Set = Year_Set,Years2Include = Years2Include,
Rotate=MapDetails_List[["Rotate"]],Cex=MapDetails_List[["Cex"]],
Legend=MapDetails_List[["Legend"]],zone=MapDetails_List[["Zone"]],
mar=c(0,0,2,0), oma=c(3.5, 3.5, 0,0), cex=1.8)
SpatialDeltaGLMM::PlotAniso_Fn(FileName=paste0(DateFile,
"Aniso.png"), Report=Report, TmbData=TmbData)
Cov_List=Summarize_Covariance(Report=Report, ParHat=Obj$env$parList(),
Data=TmbData,SD=Opt$SD, plot_cor=TRUE,
category_names=levels(Data_Geostat[,"spp"]),
plotdir=DateFile, plotTF=FieldConfig, mgp=c(2,0.5, 0),
tck=-0.02, oma=c(0,5,2,2))
PlotMap_Fn(MappingDetails, Mat, PlotDF, MapSizeRatio = c(`Width(in)` = 4,
`Height(in)` = 4), Xlim, Ylim, FileName = paste0(getwd(), "/"), Year_Set,
Rescale = FALSE, Rotate = 0, Format = "png", Res = 200, zone = NA,
Cex = 0.01, textmargin = "", add = FALSE, pch = 15,
outermargintext = c("Eastings", "Northings"), zlim = NULL, Col = NULL,
Legend = list(use = FALSE, x = c(10, 30), y = c(10, 30)), mfrow = c(1, 1),
plot_legend_fig = TRUE, land_color = "grey", ignore.na = FALSE)
SpatialDeltaGLMM::PlotResultsOnMap_Fn(plot_set = c(3),
MappingDetails = MapDetails_List[["MappingDetails"]],
Report = Report, Sdreport = Opt$SD, PlotDF = MapDetails_List[["PlotDF"]],
MapSizeRatio = MapDetails_List[["MapSizeRatio"]],
Xlim = MapDetails_List[["Xlim"]], Ylim = MapDetails_List[["Ylim"]],
FileName = DateFile, Year_Set = Year_Set, Years2Include = Years2Include,
Rotate = MapDetails_List[["Rotate"]], Cex = MapDetails_List[["Cex"]],
Legend = MapDetails_List[["Legend"]], zone = MapDetails_List[["Zone"]],
mar = c(0, 0, 2, 0), oma = c(3.5, 3.5, 0, 0), cex = 1.3,
category_names = levels(Data_Geostat[, "spp"]),land_color="grey")
Index = SpatialDeltaGLMM::PlotIndex_Fn(DirName = DateFile,
TmbData = TmbData, Sdreport = Opt[["SD"]], Year_Set = Year_Set,
Years2Include = Years2Include, strata_names = "", cex=.5,plot_legend=FALSE,
use_biascorr = TRUE,cex=0.5,category_names =c("Aequorea", "Aurelia" ,"Chrysaorea" , "Cyanea" , "Staurophora"))
pander::pandoc.table(Index$Table[, c("Category","Year", "Estimate_metric_tons", "SD_mt")])
SpatialDeltaGLMM::Plot_range_shifts(Report = Report, TmbData = TmbData, Sdreport = Opt[["SD"]],
Znames = colnames(TmbData$Z_xm), PlotDir = DateFile,
category_names = levels(Data_Geostat[, "spp"]),
Year_Set = Year_Set, mar = c(2, 2, 2,2), oma = c(1, 2, 0, 1))
Plot_Overdispersion(filename1 = paste0(DateDir, "Overdispersion"), filename2 = paste0(DateDir, "Overdispersion--panel"), Data = TmbData, ParHat = ParHat, Report = Report, ControlList1 = list(Width = 5, Height = 10, Res = 200, Units = "in"), ControlList2 = list(Width = TmbData$n_c, Height = TmbData$n_c, Res = 200, Units = "in"))
Plot_factors(Report = Report, ParHat = Obj$env$parList(), Data = TmbData, SD = Opt$SD, mapdetails_list = MapDetails_List, Year_Set = Year_Set, category_names = levels(temp[, "Sci"]), plotdir = DateFile)