Skip to content

Commit

Permalink
Expand CQL vignette with instructions on running Bayesian models with…
Browse files Browse the repository at this point in the history
… oxcAAR
  • Loading branch information
joeroe committed Sep 30, 2020
1 parent ba606ac commit 0ad2d33
Showing 1 changed file with 62 additions and 49 deletions.
111 changes: 62 additions & 49 deletions vignettes/cql.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -28,59 +28,13 @@ This vignette describes how to use this interface to generate CQL models in R.
* Execute CQL scripts in OxCal using `write_oxcal()` or the [oxcAAR package](https://CRAN.R-project.org/package=oxcAAR).


## Writing CQL scripts
## Writing CQL models

```{r setup}
library("stratigraphr")
```

## Executing CQL scripts

You can run models generated by `cql()` using the desktop or online versions of OxCal by simply copying the output into the program.
Alternatively, use `write_oxcal()` to create a .oxcal file:

```{r write-oxcal}
oxcal_cql <- cql(
cql_r_date("ABC-001", 9100, 30),
cql_r_date("ABC-002", 9200, 30),
cql_r_date("ABC-003", 9300, 30)
)
write_oxcal(oxcal_cql, "cql.oxcal")
```
```{r write-oxcal-cleanup, include=FALSE}
fs::file_delete("cql.oxcal")
```

You can also run OxCal directly through R using the [oxcAAR](https://CRAN.R-project.org/package=oxcAAR) package.
This depends on a local installation of OxCal.
If you already have one installed, you can set the path to the executable using `oxcAAR::setOxcalExecutablePath()`.
Otherwise, use `oxcAAR::quickSetupOxcal()` to download one, for example to a temporary directory:

```{r oxcaar-setup, results=FALSE, message=FALSE}
library("oxcAAR")
quickSetupOxcal(path = fs::path_temp())
```

You can then use `oxcAAR::executeOxcalScript()` to run our CQL script.
You will need to capture the path to the output file and then parse it with `oxcAAR::readOxcalOutput()` and `oxcAAR::parseOxcalOutput()`.

```{r oxcaar-execute, results=FALSE, message=FALSE}
executeOxcalScript(oxcal_cql) %>%
readOxcalOutput() %>%
parseOxcalOutput() ->
oxcal_output
```

The result is an `oxcAARCalibratedDatesList` object containing the output returned by OxCal.
Simple models (such as a list of calibrated dates) can be visualised using the built-in plotting functions in oxcAAR:

```{r oxcaar-plot, fig.show="hold"}
plot(oxcal_output)
calcurve_plot(oxcal_output)
```

## Generating CQL scripts from other types of data
## Generating CQL models from other types of data

Used in the simple way above, `stratigraphr`'s CQL interface offers little benefit over writing CQL directly.
Its real power is in combining `cql()` with other R tools to build models based on other data.
Expand Down Expand Up @@ -123,11 +77,70 @@ shub1_radiocarbon %>%
arrange(desc(phase)) %>%
summarise(cql = cql_sequence("Shubayqa 1", cql, boundaries = TRUE)) %>%
pluck("cql") %>%
cql()
cql() ->
shub1_cql
shub1_cql
```

### Stratigraphs

...

## Running CQL models

You can run models generated by `cql()` using the desktop or online versions of OxCal by simply copying the output into the program.
Alternatively, use `write_oxcal()` to create a .oxcal file:

```{r write-oxcal}
oxcal_cql <- cql(
cql_r_date("ABC-001", 9100, 30),
cql_r_date("ABC-002", 9200, 30),
cql_r_date("ABC-003", 9300, 30)
)
write_oxcal(oxcal_cql, "cql.oxcal")
```
```{r write-oxcal-cleanup, include=FALSE}
fs::file_delete("cql.oxcal")
```

You can also run OxCal directly through R using the [oxcAAR](https://CRAN.R-project.org/package=oxcAAR) package.
This depends on a local installation of OxCal.
If you already have one installed, you can set the path to the executable using `oxcAAR::setOxcalExecutablePath()`.
Otherwise, use `oxcAAR::quickSetupOxcal()` to download one, for example to a temporary directory:

```{r oxcaar-setup, results=FALSE, message=FALSE}
library("oxcAAR")
quickSetupOxcal(path = fs::path_temp())
```

You can then use `oxcAAR::executeOxcalScript()` to run the CQL script and `oxcAAR::readOxcalOutput()` to read the output back into R.

```{r oxcaar-execute, results=FALSE, message=FALSE}
executeOxcalScript(oxcal_cql) %>%
readOxcalOutput() ->
oxcal_output
```

You can parse the output with `oxcAAR::parseOxcalOutput()` and visualise it using `oxcAAR`'s built-in plotting functions:

```{r oxcaar-plot, fig.show="hold"}
oxcal_parsed <- oxcAAR::parseOxcalOutput(oxcal_output)
plot(oxcal_parsed)
calcurve_plot(oxcal_parsed)
```

The current CRAN version of oxcAAR (v. 1.0.0) does not read the posterior probabilities produced by a model with Bayesian calibration, so to work with these you need to install the latest development version (`devtools::install_github("ISAAKiel/oxcAAR")`).
With this, `oxcAAR::parseOxcalOutput()` also contains the modelled results in `$posterior_sigma_ranges` and `$posterior_probabilities`.
Again, you can quickly visualise these with the built-in plotting functions:

```{r oxcaar-bayes}
# Not run: slow
# shub1_oxcal <- executeOxcalScript(shub1_cql)
# readOxcalOutput(shub1_oxcal) %>%
# parseOxcalOutput() %>%
# plot()
```

## References

0 comments on commit 0ad2d33

Please sign in to comment.