.Rhistory

intall.packages("ggplot2")
install.packages("ggplot2")
install.packages("rtools")
seven <- 7
seven
anything
anything <- 19
anything
anything and seven
anything + seven
mean(anything + seven)
mean(anything:seven)
anything:seven
mean(c(1,2,3,4,5))
13
89
Sys.which("make")
find_rtools()
find_rtools()
load("~/data-carpentry/data-carpentry.Rproj")
156/20
156/7
##############################################################################################
#                                    LOADING LIBRARIES                                       #
##############################################################################################
library(tidyverse)
##############################################################################################
#                                    LOADING LIBRARIES                                       #
##############################################################################################
library(tidyverse)
library(rcarbon)
##############################################################################################
#                                    LOADING LIBRARIES                                       #
##############################################################################################
library(tidyverse)
install.packages("broom", type="binary")
##############################################################################################
#                                    LOADING LIBRARIES                                       #
##############################################################################################
library(tidyverse)
##############################################################################################
#                                    LOADING LIBRARIES                                       #
##############################################################################################
library(tidyverse)
##############################################################################################
#                                    LOADING LIBRARIES                                       #
##############################################################################################
library(tidyverse)
##############################################################################################
#                                    LOADING LIBRARIES                                       #
##############################################################################################
library(tidyverse)
##############################################################################################
#                                    LOADING LIBRARIES                                       #
##############################################################################################
library(tidyverse)
library(rcarbon)
library(rnaturalearth)
library(rnaturalearthhires)
library(sf)
library(sp)
##############################################################################################
#                                    PREPARING THE DATA                                      #
##############################################################################################
#-----------------------------------------Creation and management of data frame for all radiocarbon dates
c14dates <- read.csv("data_raw/c14dates.csv")                #reading the CSV file with radiocarbon dates
getwd()
dir("carib-model")
setwd("carib-model")
library(ggplot2)
library(rcarbon)
library(rnaturalearth)
library(rnaturalearthhires)
library(sf)
library(sp)
##############################################################################################
#                                    PREPARING THE DATA                                      #
##############################################################################################
#-----------------------------------------Creation and management of data frame for all radiocarbon dates
c14dates <- read.csv("data_raw/c14dates.csv")                #reading the CSV file with radiocarbon dates
c14dates$site <- factor(c14dates$site)                       #converting the column 'dates' in a factor
c14dates$country_state <- factor(c14dates$country_state)     #converting the column 'country_state' in a factor
c14dates$region <- factor(c14dates$region)                   #converting the column 'region' in a factor
c14dates$culture <- factor(c14dates$culture)                 #converting the column 'culture' in a factor
##############################################################################################
#                             CALIBRATING RADIOCARBON DATES                                  #
##############################################################################################
#--------------------------Calibrating the radiocarbon dates with normalised parameter
c14calib_normal <- calibrate(
x = c14dates$age,        #The radiocarbon age for each of the 135 radiocarbon dates
errors = c14dates$std,   #Standard deviation of each radiocarbon date
method = "mixed",        #Method for mixing the radiocarbon atmospheric curves
mixed.curves = TRUE,     #Enabling the parameter for mixing the atmospheric curves
intcal20 = intcal20(),   #Calling the IntCal20 atmospheric curve for Northern Hemisphere
shcal20 = shcal20(),     #Calling the ShCal20 atmospheric curve for Southern Hemisphere
normalised = TRUE        #Parameter for normalizing the data
)
#----------------------------Calibrating the radiocarbon dates with unnormalised parameter
c14calib_unormal <- calibrate(
x = c14dates$age,          #The radiocarbon age for each of the 135 radiocarbon dates
errors = c14dates$std,     #Standard deviation of each radiocarbon date
method = "mixed",          #Method for mixing the radiocarbon atmospheric curves
mixed.curves = TRUE,       #Enabling the argument for mixing the atmospheric curves
intcal20 = intcal20(),     #Calling the IntCal20 atmospheric curve for Northern Hemisphere
shcal20 = shcal20(),       #Calling the ShCal20 atmospheric curve for Southern Hemisphere
normalised = FALSE
)
#---------Creating a Time Range object for running the analyses
timeRange <- c(
2300,  #Starting point on 2300 BP (BC 300)
300    #Ending point on 300 BP (AD 1700)
)
###############################################################################################
#                                       #   FIRST MODEL   #                                   #
#     Summed Probability Distribution (SPD) and a Growth Model for all radiocarbon dates      #
###############################################################################################
###############################################################################################
#---------------------------  Step 1 - Creating the SPD Model  -------------------------------#
###############################################################################################
c14calib_spd <- spd(         #Generation of an object for the SPD
c14calib_unormal,          #Using the unnormalised radiocarbon calibrations
timeRange = timeRange      #Time Range between 2200 BP to 100 BP
)
###############################################################################################
#---------------------------   Step 2 - Binning the SPD Model    -----------------------------#
###############################################################################################
c14dates_bins <- binPrep(    #Generation of an object for the bins
sites = c14dates$site,     #Column for archaeological sites
ages = c14dates$age,       #Column with the radiocarbon age values
h = 100                    #Parameter for binning between an interval of 100 years
)
c14dates_binned_spd <- spd(  #Generation of the final binned SPD model
c14calib_unormal,          #Object with all the unnormalised radiocarbon dates
bins = c14dates_bins,      #Object with the bins values for binning the model
timeRange = timeRange,     #Time Range of the analysis frmo 2200 to 100 BP
)
###############################################################################################
#--------------  Step 3 - Generating a Kernel Density Estimation Model (KDE)  ----------------#
###############################################################################################
C14randates <- sampleDates(
c14calib_unormal,
bins = c14dates_bins,
nsim = 100,
verbose = FALSE
)
c14dates_ckde <- ckde(
C14randates,
timeRange = timeRange,
bw = 100
)
###############################################################################################
#--------------------- Step 4 - Generating an Exponential Growth Model ----------------------- #
###############################################################################################
c14dates_expnull <- modelTest(
c14calib_unormal,
errors = c14dates$std,
bins = c14dates_bins,
nsim = 100,
timeRange = c(2500,0),
model = "exponential",
runm = 100
)
###############################################################################################
#-------------------   Step 5 - Permutation Test of the SPD by Pottery   ---------------------#
###############################################################################################
permc14dates_culture <- permTest(
x = c14calib_unormal,
marks = c14dates$culture,
timeRange = c(2300,100),
bins = c14dates_bins,
nsim = 100,
runm = 50
)
permc14_culture_stack_spd <- stackspd(
x = c14calib_unormal,
group = c14dates$culture,
timeRange = c(2300,0),
bins = c14dates_bins,
runm = 50,
verbose = FALSE
)
###############################################################################################
#-------------------   Step 6 - Permutation Test of the SPD by region   ----------------------#
###############################################################################################
permc14dates_region <- permTest(
x = c14calib_unormal,
marks = c14dates$region,
timeRange = c(2300,100),
bins = c14dates_bins,
nsim = 100,
runm = 50
)
permc14_region_stack_spd <- stackspd(
x = c14calib_unormal,
group = c14dates$region,
timeRange = c(2300,0),
bins = c14dates_bins,
runm = 50,
verbose = FALSE
)
######################################################################################################
#-------------------------------    Step 7 - Spatial Permutation Test    ----------------------------#
######################################################################################################
allsites <- unique(data.frame(id = c14dates$site, lat = c14dates$lat, long = c14dates$long))
allsites <- st_as_sf(allsites, coords = c('long','lat'), crs = 4326)
breaks <- seq(1100,500,-100)
timerange_spatialperm <- c(1100,500)
carib_spd <- spd(
x = c14calib_unormal,
bins = c14dates_bins,
timeRange = timerange_spatialperm
)
spatial_perm <- sptest(
calDates = c14calib_unormal,
bins = c14dates_bins,
timeRange = timerange_spatialperm,
locations = allsites,
locations.id.col = 'id',
h = 100,
kernel = 'gaussian',
permute = 'locations',
nsim = 100,
breaks = breaks,
ncores = 1,
verbose = FALSE
)
world <- ne_countries(scale = 'medium', type = 'map_units', returnclass = 'sf')
xrange <- st_bbox(allsites)[c(1,3)]
yrange <- st_bbox(allsites)[c(2,4)]
######################################################################################################
#                                     #    SECOND MODEL   #                                          #
#         Generating a Summed Probability Distribution (SPD) only for Santarém pottery subset        #
######################################################################################################
c14santarem <- subset(c14dates, culture == "Santarem") #Extracting a subset of Santarém Radiocaron dates
c14calib_santarem <- calibrate(
x = c14santarem$age,          #The radiocarbon age for each of the 135 radiocarbon dates
errors = c14santarem$std,     #Standard deviation of each radiocarbon date
method = "mixed",             #Method for mixing the radiocarbon atmospheric curves
mixed.curves = TRUE,          #Enabling the parameter for mixing the atmospheric curves
intcal20 = intcal20(),        #Calling the IntCal20 atmospheric curve for Northern Hemisphere
shcal20 = shcal20(),          #Calling the ShCal20 atmospheric curve for Southern Hemisphere
normalised = FALSE
)
#Plot 01 - Summed Probability Distribution (SPD) for all radiocarbon datings
plot(
c14calib_spd,           #SPD object for all Radiocarbon dates
calendar = 'BCAD',      #Plot the SPD Model in BC/AD years
main = "Summed Probability Distribution of Koriabo, Santarém and Konduri" #Label
)
plot(                     #Plotting the smoothing line with a bin value of 200 years
c14calib_spd,           #SPD object for all radiocarbon dates
calendar = 'BCAD',      #Calendar in BC/AD years
runm = 200,             #Smoothed distribution with a binning of 200 years
add = TRUE,             #Add the line
type = "simple",        #Simple line
col = "red",            #Line colour = Red
lwd = 2,                #Line Thickness = 2
lty = 2
)
#Plot 02 - Binned Summed Probability Distribution (SPD) for all radiocarbon datings
plot(
c14dates_binned_spd,    #SPD object for all Radiocarbon dates
calendar = 'BCAD',      #Plot the SPD Model in BC/AD years
main = "Summed Probability Distribution of Koriabo, Santarém and Konduri" # Plot Title
)
plot(                     #Plotting the smoothing line with a bin value of 200 years
c14dates_binned_spd,    #SPD object for all radiocarbon dates
calendar = 'BCAD',      #Calendar in BC/AD years
runm = 200,             #Smoothed distribution with a binning of 200 years
add = TRUE,             #Add the line
type = "simple",        #Simple line
col = "red",            #Line colour = Red
lwd = 2,                #Line Thickness = 2
lty = 2
)
#Plot XX - Kernel Density Estimation (KDE) for all dates
plot(
c14dates_ckde,
type = 'multiline',
main = "Kernel Density Estimation (KDE) for Koriabo, Konduri and Santarém"
)
#Plot XX - Exponential Growth Model for all dates
plot(c14dates_expnull,
calendar = 'BCAD',
main = "Exponential Growth Model (Koriabo, Konduri, Santarém)"
)
#Plot XX - Permutation Test for Koriabo, Santarém and Konduri cultures
par(mfrow = c(1,1))
plot(permc14dates_culture, focalm = 1, calendar = 'BCAD', main = "Koriabo")
plot(permc14dates_culture, focalm = 2, calendar = 'BCAD', main = "Santarém")
plot(permc14dates_culture, focalm = 3, calendar = 'BCAD', main = "Konduri")
par(mfrow=c(2,2))
plot(permc14_culture_stack_spd, type = 'stacked', calendar = 'BCAD')
plot(permc14_culture_stack_spd, type = 'lines',  calendar = 'BCAD')
plot(permc14_culture_stack_spd, type = 'multipanel',  calendar = 'BCAD')
plot(permc14_culture_stack_spd, type = 'proportion',  calendar = 'BCAD')
#Plot XX - Permutation Test for 5 different Regions
par(mfrow = c(2,3))
plot(permc14dates_region, focalm = 1, calendar = 'BCAD', main = "Western Guiana")
plot(permc14dates_region, focalm = 2, calendar = 'BCAD', main = "Guiana Coast")
plot(permc14dates_region, focalm = 3, calendar = 'BCAD', main = "Amapa-Estuary")
plot(permc14dates_region, focalm = 4, calendar = 'BCAD', main = "Lower Amazon")
plot(permc14dates_region, focalm = 5, calendar = 'BCAD', main = "Xingu-Iriri Rivers")
par(mfrow=c(2,2))
plot(permc14_region_stack_spd, type = 'stacked', calendar = 'BCAD')
plot(permc14_region_stack_spd, type = 'lines',  calendar = 'BCAD')
plot(permc14_region_stack_spd, type = 'multipanel',  calendar = 'BCAD')
plot(permc14_region_stack_spd, type = 'proportion',  calendar = 'BCAD')
#Plot XX - Block of Growth Rates
plot(spd2rc(carib_spd, breaks = breaks))
#Plot XX - Spatial Permutation Test
par(mar=c(2,2,2,2), mfrow=c(2,1))
#Plot1
plot(world$geometry,
col = "antiquewhite3",
border = "gray20",
xlim = xrange,
ylim = yrange,
main = "Growth Rate "
)
plot(spatial_perm,
index = 2,
option = "raw",
add = TRUE,
breakRange = c(-0.005,0.005),
breakLength = 7,
rd = 5,
legend = TRUE,
legSize = 0.7,
baseSize = 1.7,
location = "topright")
#Plot2
plot(world$geometry,
col = "antiquewhite3",
border = "gray20",
xlim = xrange,
ylim = yrange,
main = "Growth Rate "
)
plot(spatial_perm,
index = 3,
option = "raw",
add = TRUE,
breakRange = c(-0.005,0.005),
breakLength = 7,
rd = 5,
legend = TRUE,
legSize = 0.7,
baseSize = 1.7,
location = "topright")
#Plot3
plot(world$geometry,
col = "antiquewhite3",
border = "gray20",
xlim = xrange,
ylim = yrange,
main = "Growth Rate "
)
plot(spatial_perm,
index = 3,
option = "raw",
add = TRUE,
breakRange = c(-0.005,0.005),
breakLength = 7,
rd = 5,
legend = TRUE,
legSize = 0.7,
baseSize = 1.7,
location = "topright")
#Plot4
plot(world$geometry,
col = "antiquewhite3",
border = "gray20",
xlim = xrange,
ylim = yrange,
main = "Growth Rate "
)
plot(spatial_perm,
index = 4,
option = "raw",
add = TRUE,
breakRange = c(-0.005,0.005),
breakLength = 7,
rd = 5,
legend = TRUE,
legSize = 0.7,
baseSize = 1.7,
location = "topright")
#Plot5
plot(world$geometry,
col = "antiquewhite3",
border = "gray20",
xlim = xrange,
ylim = yrange,
main = "Growth Rate "
)
plot(spatial_perm,
index = 5,
option = "raw",
add = TRUE,
breakRange = c(-0.005,0.005),
breakLength = 7,
rd = 5,
legend = TRUE,
legSize = 0.7,
baseSize = 1.7,
location = "topright")
#Plot6
plot(world$geometry,
col = "antiquewhite3",
border = "gray20",
xlim = xrange,
ylim = yrange,
main = "Growth Rate "
)
plot(spatial_perm,
index = 6,
option = "raw",
add = TRUE,
breakRange = c(-0.005,0.005),
breakLength = 7,
rd = 5,
legend = TRUE,
legSize = 0.7,
baseSize = 1.7,
location = "topright")
install.packages("tidyverse")
library(tidyverse)