Jeffrey Stevens

In the last post, we used the qualtRics

 1

package to import survey data directly from Qualtrics accounts.I often use Qualtrics for survey data collected via Amazon’s Mechanical Turk, which is an automated system that connects online respondents to paid surveys.However, sometimes those data are not super useful because the quality of the data is poor.In the anonymous, large-scale market of online survey respondents, the financial

Working with Qualtrics Data - Part 2: Excluding Data

Erik Bülow

<figure class="pull-right">
<a href="https://docs.ropensci.org/coder/"><img src="hex.png"
width="150" alt="Hex sticker for the coder package" /></a>
</figure>

A new R-package, [coder](https://docs.ropensci.org/coder/), has been
[developed](https://github.com/ropensci/coder),
[peer-reviewed](https://github.com/ropensci/software-review/issues/381)
by rOpenSci, accepted by
[CRAN](https://cran.r-project.org/package=coder), and published in a
paper by the Journal of Open Source Software
[(JOSS)](https://joss.theoj.org/papers/10.21105/joss.02916). In this
blog post, I will explain why this package might be useful for
(epidemiological/medical/health care related) research.

## Clinical mess

Once upon a time, in countries not far from ours, there were MDs and
nurses making up funny names for any diseases they encountered. What Dr.
A called X would be recognized only as Y by his dear colleague Dr. B. X,
however, was also a name used by Dr. C, but then for a completely
different condition. It was a mess!

Those times are long gone due to the fascinating story of standardized
medical nomenclatures and coding. It all started in the first half of
the 19th century and has been well documented in a [free on-line
book](https://www.cdc.gov/nchs/data/misc/classification_diseases2011.pdf)
by Iwao Moryama, Ruth Loy and Alastair Robb-Smith.

Today, when visiting the hospital, our (somatic) conditions are most
likely recorded in a medical register by some adaptation of the 10th
revision of the International Statistical Classification of Diseases and
Related Health Problems (ICD-10). This is a global standard for medical
coding administrated by the World Health Organization (WHO). Psychiatric
conditions might as well be recorded by the Diagnostic and Statistical
Manual of Mental Disorders (DSM). Prescribed medical substances are
reported by the Anatomical Therapeutic Chemical (ATC) classification
system. Surgical procedures performed in the Nordic countries are coded
by the Nordic Medico-Statistical Committee's (NOMESCO's) Classification
of Surgical Procedures (NCSP). All combined, there are several hundred
thousands of different codes quantifying every possible interaction with
health care. Although the mess has significantly decreased, the world of
medical coding is still rather messy! Not only are there different
versions of ICD (currently with ICD-11 in the pipe). There are also
yearly revisions of the currently used version, as well as clinical
versions, oncological versions, and several national implementations.
All this might be well motivated from a clinical and/or administrative
view-point, but it does make things difficult for researchers trying to
pool data from different sources, countries and periods.

## Comorbidity

A measure of patient comorbidity is used to summarize the total
pre-existing disease burden prior to other diagnosis or events of
interest. Such data are often used to control for case-mix confounding
in statistical models with various exposures and outcomes. There is no
obvious definition of what to consider a comorbid condition, but several
suggestions exist. [Charlson et
al.](https://pubmed.ncbi.nlm.nih.gov/3558716/), [Elixhauser et
al.](https://pubmed.ncbi.nlm.nih.gov/9431328/), [Sloan et
al.](https://pubmed.ncbi.nlm.nih.gov/12773842/) and others, have
combined individual medical conditions into more holistic
categories/diseases. The number of categories might then be summarized
by an index value (a weighed sum), as a single comorbidity-measure for
each patient. There are several versions of each index, however. To use
those indices in practice is also cumbersome since many data sources
(large administrative registers) are big - bigger than what would fit in
the random access memory (RAM) of a standard computer.

## R packages

There are some excellent R-packages to help with both the Charlson and
Elixhauser comorbidity indices.
[comorbidity](https://ellessenne.github.io/comorbidity/) by Alessandro
Gasparini and [icd](https://jackwasey.github.io/icd/) by Jack O. Wasey
and Michael Lang are very well documented with efficient
implementations. To keep up with all different and newly proposed coding
versions seems like a daunting, if not impossible task, however.
Therefore, the coder package offers an alternative, more flexible,
implementation. It is not hard-wired to Charlson or Elixhauser, although
capabilities for those indices are included as well.
[data.table](https://rdatatable.gitlab.io/data.table/) is used
internally for computational efficiency, while the application
programming interface (API) is more inspired by
[tidyverse](https://www.tidyverse.org/).

## Regular expressions

In coder, all relevant codes (from ICD, ATC, NCSP or whatever) are
represented in a compact way by their corresponding regular expressions.
Peptic ulcer disease for example is recognized as an important
comorbidity by the Elixhauser classification. It might be coded by
ICD-10 as "K25.x--K28.x" (see [table 2
here](https://journals.lww.com/lww-medicalcare/Abstract/2005/11000/Coding_Algorithms_for_Defining_Comorbidities_in.10.aspx)).
This format is well understood by humans; "-" indicates an interval and
"x" acts like a wild card for any additional alphanumerical characters.
This might be translated into a code list (ignoring the dots): K257,
K259, K267, K269, K277, K279, K287, K289. But why those codes only? Why
not K250 or K289A? Well, not all codes are used in practice. This list
is just an example corresponding to codes defined in the clinical
implementation of ICD-10 (ICD-10-CM) used in 2020 in the USA. The list
might look different with data from the Swedish version of ICD-10
(ICD-10-SE) recorded in 2021. Hence, the use of regular expressions will
relax the need of a complete and version-specific code list.A regular
expression for peptic ulcer disease might be `^K2[5-8][79]$` where `^`
(`$`) marks the beginning (end) of the character string, where
alternative digits are written within brackets, and where the generic
wild card `x` has been limited to either 7 or 9 as used in clinical
practice (see `?regex` in R for details).

To use regular expressions is also more computationally efficient
compared to a traditional approach where each individual code might be
compared by exact matches to any "word" (cell) in a long dictionary
(vector) of individual codes. Hence, it is easy to realize that the code
`K257` can not be found in a dictionary limited to words with initials
`L` to `Z`. Thus, you do not need to compare the code `K257` to each and
every entry in that list to see if you get a match or not.Such
computational efficiency is crucial if working with large populations
(as from a national patient register spanning several decades et
cetera).

By default, there are 413 of those regular expressions included in the
coder package. They cover all conditions recognized by Charlson,
Elixhauser, RxRisk V, the comorbidity-polypharmacy score (CPS), as well
as for some diagnose specific adverse events classifications after hip
or knee replacement. Charlson, Elixhauser and Rx Risk V have 6, 5 and 3
alternative regular expressions for each condition. Peptic ulcer disease
for example can also be recognized from 1) `53([1-4][79]0)|V1271`; 2)
`53([1-4]([4-69]1|7))`; or 3) `53[1-4][79]` based on three different
versions of ICD-9-CM. Additional versions (think of the upcoming ICD-11)
are relatively easy to implement by using a dedicated `classcodes`
object (an enhanced data frame/tibble with alternative regular
expressions in each column).

## Interpretation

Although `classcodes` objects are relatively easy to implement, the
interpretation of a regular expression might be hard to grasp
intuitively.Or what about hypertension according to [Pratt's
version](https://bmjopen.bmj.com/content/8/4/e021122) of the Rx Risk V
classification based on ATC codes:
`C0(2(A(B[0-9]{2}|C0[0-5])|DB(0[2-9]|[1-9][0-9]))|3((A[A-Z][0-9]{2}|BA(0[0-9]|1[01]))|DB(01|99)|C([AB][0-9]{2}|C0[01])|EA01)|9(BA0[2-9]|DA0[2-8]|C[A-X][0-9]{2}))`?Fortunately,
there are several functions in the coder package to aid such
interpretation:

- `visualize()` will open the default web browser to show a structural
 diagram based on the [JavaScript Regular Expression
 Visualizer](https://jex.im/regulex).
- `summary()` relies on the
 [decoder](https://cran.r-project.org/package=decoder) package to
 decode the meaning of a regular expression in the content of a
 specified code version. For example the ICD-10 specification of a
 dislocation after total hip arthroplasty (`^(M24(3|4F?)|S730|T933)$`)
 would expand to 5 individual codes according to ICD-10-SE: M243, M244,
 M244F, S730, T933. None of those would be recognized by ICD-10-CM,
 however. A longer list of 146 codes (M2430, M24311, M24312, M24319,
 M24321, ...) would be recognized for `^(M24(3|4F?)|S730|T933)`
 (without the `$` suffix).
- `codebook()` takes this a step further. It not only lists the relevant
 codes but also translates them into plain text.M2430 and M24311 from
 ICD-10-CM for example would be described as "Pathological dislocation
 of unspecified joint, not elsewhere classified" and "Pathological
 dislocation of right shoulder, not elsewhere classified". A code such
 as C430 would be recognized by both ICD-10-CM and ICD-10-SE, with the
 same meaning but in different languages: "Malignant melanoma of lip"
 and "Malignt melanom på läpp".

## Use case

Assume we have a list of patients who received surgery at specified
dates:

::: {.highlight}
``` {.chroma tabindex="0"}
library(coder)ex_people # Random names and dates included in the package
```
:::

``` {tabindex="0"}
# A tibble: 100 x 2 name surgery <chr> <date> 1 Chen, Trevor 2021-01-19 2 Graves, Acineth 2020-10-11 3 Trujillo, Yanelly 2020-09-28 4 Simpson, Kenneth 2020-12-31 5 Chin, Nelson 2020-12-14 6 Le, Christina 2020-07-18 7 Kang, Xuan 2020-10-20 8 Shuemaker, Lauren 2020-07-19 9 Boucher, Teresa 2020-12-2510 Le, Soraiya 2020-11-29# ... with 90 more rows
```

Those patients (among others) were also recorded in a national patient
register with date of hospital admissions and corresponding ICD-10
codes. Each patient can have multiple hospital visits over the years.
Also, multiple diagnosis might have been recorded at each occasion.

::: {.highlight}
``` {.chroma tabindex="0"}
ex_icd10
```
:::

``` {tabindex="0"}
# A tibble: 2,376 x 4 name admission icd10 hdia <chr> <date> <chr> <lgl> 1 Tran, Kenneth 2020-08-02 S134A FALSE 2 Tran, Kenneth 2021-01-16 W3319 FALSE 3 Tran, Kenneth 2020-12-26 Y0262 TRUE 4 Tran, Kenneth 2020-11-18 X0488 FALSE 5 Sommerville, Dominic 2021-01-07 V8104 FALSE 6 Sommerville, Dominic 2020-08-18 B853 FALSE 7 Sommerville, Dominic 2021-01-02 Q174 FALSE 8 Sommerville, Dominic 2020-08-23 A227 FALSE 9 Sommerville, Dominic 2020-12-28 H702 FALSE10 Sommerville, Dominic 2020-04-21 X6051 TRUE # ... with 2,366 more rows
```

Using those two data sets, as well as a classification scheme (a
`classcodes` object, here `charlson`), we can easily identify all
Charlson comorbidities for each patient:

::: {.highlight}
``` {.chroma tabindex="0"}
ch <- categorize( ex_people, # patients of interest codedata = ex_icd10, # Medical codes from national patient register cc = "charlson", # Calculate Charlson comorbidity id = "name", code = "icd10" # Identify column names )
```
:::

``` {tabindex="0"}
Classification based on: icd10
```

Now, how many of those patients were diagnosed with malignancy?

::: {.highlight}
``` {.chroma tabindex="0"}
sum(ch$malignancy)
```
:::

``` {tabindex="0"}
[1] 5
```

What is the distribution of the combined Charlson comorbidity index for
each patient? (This would summarize the general health status of the
cohort.)

::: {.highlight}
``` {.chroma tabindex="0"}
barplot(table(ch$charlson))
```
:::

<figure>
<img src="unnamed-chunk-5-1.png"
alt="Bar plot showing that most patients have 0 comorbidity." />
<figcaption>Distribution of the Charlson comorbidity index. Patients
in this cohort are healthy. Most have zero comorbidities and only a few
have several.</figcaption>
</figure>

## Further usage

There are many more examples provided at the [coder web
page](https://docs.ropensci.org/coder/). It is easy to constrain the
time frame to comorbidities recorded during only one year prior to
surgery, or to 90 days for adverse events thereafter, or to use
different versions of the regular expressions etc. Core functionality is
provided by the `categorize()` function used above. Alternatively, more
fine grained control is provided through a more explicit chain of
commands: `codify() %>% classify() %>% index()`.

So far, the package has been used in several research projects (cited in
the corresponding [JOSS
article](https://joss.theoj.org/papers/10.21105/joss.02916)).

Contributions are most welcome and I would be very happy to [hear
about](https://github.com/ropensci/coder/issues) possible additional use
cases. And please note that there are no hard-coded limitations to
patients and/or medical data. The package can be used with all types of
codified data!

## Acknowledgement

An early version of the coder package was significantly improved after
generous contributions made by [Emily C.
Zabore](https://github.com/zabore) and [David
Robinson](/author/david-robinson/) during the [the peer-review
process](https://github.com/ropensci/software-review/issues/381) for
rOpenSci. I therefore wish to express my deepest gratitude to both of
them, as well as to the rOpenSci community at large!

A new R-package, coder, has been developed, peer-reviewed by rOpenSci, accepted by CRAN, and published in a paper by the Journal of Open Source Software (JOSS). In this blog post, I will explain why this package might be useful for (epidemiological/medical/health care related) research. Clinical mess   Once upon a time, in countries not far from ours, there were MDs and nurses making up funny names for any diseases they encountered.

coder Makes Medical Coding less Messy

Lluís Revilla Sancho

In this post I will explain the history behind
[BaseSet](https://cran.r-project.org/package=BaseSet "BaseSet on CRAN")
then a brief introduction to sets, followed by showing what you can do
with [BaseSet](https://docs.ropensci.org/BaseSet/).

## Brief BaseSet history

I study diseases to try to find what causes them at a research institute
associated with an hospital.Thanks to recent technological advances we
can analyze many things from a single patient's sample.Having so much
information available can be overwhelming, making it difficult to find
the causes of diseases (but it is much better than not having enough
information!).

In order to make it easier to focus and aggregate knowledge, scientists
have been storing what we know of each gene in several ways.As more
information was gathered, efforts to group these genes into those to
focus on or discard as causes of different processes of interest started
to emerge^[1](#fn:1){.footnote-ref role="doc-noteref"}^.Scientists
categorize genes according to where they are, what do they do, how they
function, and with which other genes they interact.There are also many
different sub-groups within these categories.

This has lead to several ways of storing this information, and depending
on which database and which method you use you end up with different
groups.This in turn means that we aren't so sure of how to group genes
and as new experiments are carried out these groupings might change.

I had previously created a package
([BioCor](https://bioconductor.org/packages/BioCor/ "BioCor: Comparing sets of genes"))
to compare these gene sets which returns a numerical value of how
similar two groups of genes are.But I found out that depending on which
database I used I got very different results.So I wanted an easy way to
know in which groups a gene was present, but I also wanted this
information unified from several data sources to account for this
uncertainty.

I started [experimenting](https://github.com/llrs/GSEAdv/ "GSEAdv") by
extending a class present on
[GSEABase](https://bioconductor.org/packages/GSEABase/ "Bioconductor package")
that I had used previously.I also found out that there were some efforts
by other people through the
[Bioconductor](https://bioconductor.org "Project for high-throughput genomic data")
slack channel.So I got in touch with them and we shared some ideas on
how to do this.

To explore different ways to solve the problem we ended up developing
three different packages that aim to do similar things but take
different approaches.At the moment BaseSet published on CRAN and
[BiocSet](https://bioconductor.org/packages/BiocSet/ "BiocSet displays different biological sets in a triple tibble format")
is published on
Bioconductor.[unisets](https://github.com/kevinrue/unisets) is available
on GitHub although development has stalled due to time constraints.These
packages provide methods to work with sets, the mathematical name for
groups of elements.

## Sets

Sets are not only used for tennis games.In mathematics if you have a
group of multiple elements such as `c("apple", "bannana", "pineaple")`
it is called a set.A set is a group of elements, usually just unique
elements are used (no `c("apple", "bannana", "pineaple", "apple")`).

To learn and understand the world we group elements together.Usually it
is clear what should be inside a group, like the elements above, all of
them are fruits.Sometimes it is not so clear if an element belongs to a
certain group: [Are tomatoes
fruits?](https://en.wikipedia.org/wiki/Tomato#Fruit_versus_vegetable "Wikipedia page about tomatoes and its classification")

An element can also belong to several groups.Basketball is a team sport,
but is also played with a ball.So basketball can be in both groups, as
we see on
[Wikipedia](https://en.wikipedia.org/wiki/Basketball "Wikipedia: Basketball"):

> [Categories](https://en.wikipedia.org/wiki/Help:Category "Help:Category"):
> [Basketball](https://en.wikipedia.org/wiki/Category:Basketball "Category:Basketball")
> \| [Sports originating in the United
> States](https://en.wikipedia.org/wiki/Category:Sports_originating_in_the_United_States "Category:Sports originating in the United States")
> \| [Team
> sports](https://en.wikipedia.org/wiki/Category:Team_sports "Category:Team sports")
> \| [Summer Olympic
> sports](https://en.wikipedia.org/wiki/Category:Summer_Olympic_sports "Category:Summer Olympic sports")
> \| [Ball
> games](https://en.wikipedia.org/wiki/Category:Ball_games "Category:Ball games")
> \| [Games and sports introduced in
> 1891](https://en.wikipedia.org/wiki/Category:Games_and_sports_introduced_in_1891 "Category:Games and sports introduced in 1891")

Some say that tomatoes should be both considered as fruits and as
vegetables, so they could also be in both sets.We can do this for any
element we can think of, pencils, mathematical operations, numbers,
words, genes...

When we have more than one set we might be interested in which elements
are in two groups, which fruits are also vegetables, or which fruits are
not vegetables.These operations are called set operations.Every time you
have used `%in%` or `merge` in your code you have been using a set
operation.There are many set operations, the most commonly used are
"intersection", for which R has the `intersect` function, and "union",
which is done with `c` on vectors.The equivalents on `data.frame`s are
done with `merge` and its arguments (if you are familiar with the
`*_join` verbs from dplyr, the many arguments of `merge` are one of the
reasons why there are so many `*_join` functions).

As you see R already provides some functions to work with sets.So why
write a new package for sets?Let's see it with some code:

## Using BaseSet

First load the package:

::: {.highlight}
``` {.chroma tabindex="0"}
library("BaseSet")packageVersion("BaseSet")[1] '0.0.14'
```
:::

In order to work with sets the package defines its own class named
`TidySet`.We can transform our data with the `tidySet` function:

::: {.highlight}
``` {.chroma tabindex="0"}
sets <- list(fruits = c("apple", "bannana", "pineaple"))TS <- tidySet(sets)TS elements sets fuzzy1 apple fruits 12 bannana fruits 13 pineaple fruits 1
```
:::

We can explore some properties of this set:

::: {.highlight}
``` {.chroma tabindex="0"}
nElements(TS)[1] 3
```
:::

::: {.highlight}
``` {.chroma tabindex="0"}
nSets(TS)[1] 1
```
:::

then add more fruits

::: {.highlight}
``` {.chroma tabindex="0"}
TS1 <- add_elements(TS, "orange")nElements(TS1)[1] 4
```
:::

::: {.highlight}
``` {.chroma tabindex="0"}
TS1 elements sets fuzzy1 apple fruits 12 bannana fruits 13 pineaple fruits 1
```
:::

We can see it has one more element but it isn't displayed, that's is
because we didn't record to what set this element belongs to:

::: {.highlight}
``` {.chroma tabindex="0"}
relation <- data.frame(sets = "fruits", elements = "orange")TS1 <- add_relation(TS, relation)TS1 elements sets fuzzy1 apple fruits 12 bannana fruits 13 pineaple fruits 14 orange fruits 1
```
:::

As you might have noticed there is a third column named fuzzy that so
far is always 1.This column contains a value between 0 and 1 and
represents the membership of an element to that set.A value of 1
indicates full membership and 0 that there isn't any relationship.Values
between encode the uncertainty.We will demonstrate this later, but first
let's create a new set to do some operations with it:

::: {.highlight}
``` {.chroma tabindex="0"}
relation <- data.frame(sets = c("vegetables", "fruits"), elements = c("tomatoes", "tomatoes"))TS2 <- add_relation(TS1, relation)TS2 elements sets fuzzy1 apple fruits 12 bannana fruits 13 pineaple fruits 14 orange fruits 15 tomatoes vegetables 16 tomatoes fruits 1
```
:::

Now that we have two sets we can perform more operations:

::: {.highlight}
``` {.chroma tabindex="0"}
intersection(TS2, c("vegetables", "fruits")) elements sets fuzzy1 tomatoes vegetables∩fruits 1
```
:::

::: {.highlight}
``` {.chroma tabindex="0"}
union(TS2, c("vegetables", "fruits")) elements sets fuzzy1 apple vegetables∪fruits 12 bannana vegetables∪fruits 13 pineaple vegetables∪fruits 14 orange vegetables∪fruits 15 tomatoes vegetables∪fruits 1
```
:::

Note that the names of the resulting sets have been automatically set up
using the [appropriate
symbols](https://www.mathsisfun.com/sets/symbols.html "Link to a website with some set operations symbols").

So far we have set up how to count and merge different sets but one of
the goals of BaseSet is to also encode the uncertainty when we assign an
element to a set.

To account for uncertainty mathematicians have come up with the idea of
fuzzy sets.Fuzzy sets are sets but instead of knowing for sure that an
element belongs to the set, there is a value, the fuzzy value or
membership, that measures the uncertainty of an element belonging to a
set.

A typical example involves ratings.Imagine classification methods that,
given a picture, classify the probability of each object being in the
picture.It might return probabilities such as:

::: {.highlight}
``` {.chroma tabindex="0"}
classification <- matrix(c(0.1, 0.2, 0.6, 0.85, 0.6, 0.5, 0.25, 0.16), nrow = 2, ncol = 4, byrow = TRUE, dimnames = list(c("image1", "image2"), c("Cat", "Dog", "Duck", "Fish")))classification # Algorithm classification Cat Dog Duck Fishimage1 0.1 0.2 0.60 0.85image2 0.6 0.5 0.25 0.16
```
:::

::: {.highlight}
``` {.chroma tabindex="0"}
TS <- tidySet(classification) # Transform dataTS elements sets fuzzy1 image1 Cat 0.102 image2 Cat 0.603 image1 Dog 0.204 image2 Dog 0.505 image1 Duck 0.606 image2 Duck 0.257 image1 Fish 0.858 image2 Fish 0.16
```
:::

We can calculate how probable it is that the image is classified to just
one object.

::: {.highlight}
``` {.chroma tabindex="0"}
element_size(TS) elements size probability1 image1 0 0.04322 image1 1 0.32523 image1 2 0.48024 image1 3 0.14125 image1 4 0.01026 image2 0 0.12607 image2 1 0.38108 image2 2 0.36209 image2 3 0.119010 image2 4 0.0120
```
:::

Given this memberships there is just a 0.3252 probability that the
machine would classify the first image as just one thing, so it is more
probable that the first image will be
misclassified^[2](#fn:2){.footnote-ref role="doc-noteref"}^.

For the second image it is more probable to contain just one object than
two, but there is still high uncertainty as the probability to classify
to none of the 4 objects is 0.1260 ^[3](#fn:3){.footnote-ref
role="doc-noteref"}^.

This was not evident from the "output" of the classification but is
equally important to know in many situations.

## Finishing

I [submitted](https://github.com/ropensci/software-review/issues/359)
this package to [rOpenSci](https://ropensci.org "rOpenSci") peer-review
because I wanted more feedback on these fuzzy sets methods and wasn't
sure of the mix I made of S4 methods and S3 methods.

The topic is not a common one so although I submitted the package on
January for review it started on June (after 2 months of hiatus due to
COVID).This long waiting time occurred because it was hard for the
editor (thanks [Anna Krystalli](/author/anna-krystalli)!) to find two
reviewers with all that was going on.

Once two reviewers were found the feedback started coming (Many thanks
[Zebulun Arendsee](https://github.com/arendsee) and [Jennifer
Chang](https://github.com/j23414)!).The reviewers pointed out
inconsistencies on the order of the output, suggested more methods,
pointed out mismatches between the documentation and the code.

Most importantly much of the comments and exchange was around some
methods and explanations for fuzzy sets which were not
clear^[4](#fn:4){.footnote-ref role="doc-noteref"}^.

I made my best efforts to improve the package and hope you'll find it
useful in your projects.[Try it
out!](https://docs.ropensci.org/BaseSet/)

::: {.footnotes role="doc-endnotes"}

------------------------------------------------------------------------

1. ::: {#fn:1}
 See [this
 article](https://doi.org/10.1371/journal.pcbi.1002375 "Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges")
 reviewing the different methods
 used. [↩︎](#fnref:1){.footnote-backref role="doc-backlink"}
 :::

2. ::: {#fn:2}
 We add up the probability of being one feature (and not the others,
 i.e.: `(1- x)`) for all the elements:
 `0.1*(1-0.2)*(1-0.6)*(1-0.85) + (1-0.1)*0.2*(1-0.6)*(1-0.85) + (1-0.1)*(1-0.2)*0.6*(1-0.85) + (1-0.1)*(1-0.2)*(1-0.6)*0.85` [↩︎](#fnref:2){.footnote-backref
 role="doc-backlink"}
 :::

3. ::: {#fn:3}
 Similar to the other calculation, this probability is
 `(1-0.6)*(1-0.5)*(1-0.25)*(1-0.16)`. [↩︎](#fnref:3){.footnote-backref
 role="doc-backlink"}
 :::

4. ::: {#fn:4}
 If there is still something unclear; let me know!Open an
 [issue](https://github.com/ropensci/BaseSet/issues/new). [↩︎](#fnref:4){.footnote-backref
 role="doc-backlink"}
 :::
:::

In this post I will explain the history behind BaseSet then a brief introduction to sets, followed by showing what you can do with BaseSet.

Introducing BaseSet for mathematical sets

Laura Graham

The [grainchanger](https://docs.ropensci.org/grainchanger/) package
provides functionality for data aggregation to a coarser resolution via
moving-window or direct methods.

### Why do we need new methods for data aggregation?

As landscape ecologists and macroecologists, we often need to aggregate
data in order to harmonise datasets. In doing so, we often lose a lot of
information about the spatial structure and environmental heterogeneity
of data measured at finer resolution.

The issues around scale disconnects are both conceptual and practical:

#### Conceptual issues

- The scale at which drivers affect responses varies

- As does the resolution at which broad-scale patterns of ecological
 processes emerge

#### Technical issues

- The scale we measure this stuff at varies too

- Ecological response data (e.g. species richness from species' atlas
 data) are often available at coarse resolutions, whereas with
 increasing availability of high-resolution remote sensing data,
 environmental data are available at much finer resolutions

So how do we aggregate this fine-resolution data without losing
information about spatial structure? We argue that using a moving window
to calculate heterogeneity before aggregating retains information about
spatial structure at an appropriate scale. We can then use this
aggregated measure to (1) identify the appropriate scale-of-effect when
we only have coarse-grain response data; and (2) investigate the effect
of cross-scale interactions on ecological phenomena. See [our paper in
Methods in Ecology &
Evolution](https://doi.org/10.1111/2041-210X.13177)^[1](#fn:1){.footnote-ref
role="doc-noteref"}^ for simulations and a case study which show this.

We developed this method and R package in order to overcome some of
these issues.

### How does the method work?

The moving-window data aggregation (MWDA) method smooths an input raster
using a specified function within a moving window of a specified size
and shape prior to aggregation.

<figure>
<img src="/img/blog-images/2019-07-16-grainchanger/schematic.jpg"
width="960"
alt="Schematic of the moving-window data aggregation approach" />
</figure>

The above is a graphical representation of the MWDA method. In
calculating the MWDA measure, three aspects of scale are considered.
Predictor grain is the characteristic spatial scale of a predictor
variable, that is, the resolution of the environmental data;
scale‐of‐effect determines the appropriate scale of the relationship
between predictor and response, for example, an ecological
neighbourhood; response grain is the grain of the unit into which you
are predicting, that is, the resolution of the response variable
(represented by the black lines). Note that the colour scale is
unitless. Yellow cells represent 'high' values and dark blue cells 'low'
values. Panel 1 shows a close up of one of the response grain cells in
panel 2, whereas panel 2 shows all response grain cells for the study
region. Panel 3 shows the study region after aggregation.

The package rests around three key functions:

- `winmove`
- `winmove_agg`
- `nomove_agg`

In `winmove` the user can smooth a raster using a window of a specified
shape and size and a specified function (either those that are built in,
or a user defined function).

In `winmove_agg`, the user specifies an aggregation function in addition
to the above.

In `nomove_agg`, the user specifies only the aggregation function.

Below is an example of the three functions applied to a continuous
landscape with the moving window function as `var` within a square
window of 5 units. For `winmove_agg` the aggregation function is `mean`,
for `nomove_agg` it is `var`.

::: {.highlight}
``` {.chroma tabindex="0"}
library(grainchanger)library(ggplot2)library(landscapetools)library(patchwork)# continuous landscapep1 <- show_landscape(cont_ls, legend_title = "", xlab = "", ylab = "") + ggtitle("Landscape") + theme_bw() + coord_equal()# moving-window smoothing using variancewin_smooth <- winmove(fine_dat = cont_ls, d = 5, type = "rectangle", win_fun = var)p2 <- show_landscape(win_smooth, legend_title = "", xlab = "", ylab = "") + ggtitle("winmove") + theme_bw() + coord_equal()# moving-window aggregation using varianceg_sf$mwda <- winmove_agg(coarse_dat = g_sf, fine_dat = cont_ls, d = 5, type = "rectangle", win_fun = var, agg_fun = mean, quiet = TRUE)p3 <- ggplot(g_sf) + geom_sf(aes(fill = mwda)) + scale_fill_viridis_c(name = "") + theme_bw() + ggtitle("winmove_agg")# standard aggregation using varianceg_sf$dda <- nomove_agg(coarse_dat = g_sf, fine_dat = cont_ls, agg_fun = var, quiet = TRUE)p4 <- ggplot(g_sf) + geom_sf(aes(fill = dda)) + scale_fill_viridis_c(name = "") + theme_bw() + ggtitle("nomove_agg")(p1 + p2) / (p3 + p4)
```
:::

![](img/blog-images/2019-07-16-grainchanger/example1-1.png)

See the [package
website](https://docs.ropensci.org/grainchanger/articles/functions.html)
for a description of the built-in functions.

### Why not just use the `raster` package?

The key functionality of `grainchanger` is moving window smoothing, and
aggregation to a coarser resolution. "But aren't these available in the
`raster` package?", I hear you ask.

While `raster::focal` is super fast for calculating using `sum`, on
large rasters it can become computationally cumbersome with other
functions. However, using the weights matrix and the `sum` function, it
is possible to calculate many other functions. Our package provides a
wrapper to do so.

As a side note, it was a [helpful post from whuber on
StackExchange](https://gis.stackexchange.com/questions/151962/calculating-shannons-diversity-using-moving-window-in-r)
which led me to realise this.

Additionally, the aggregation functions in the `grainchanger` package
have been written so that (a) the user can aggregate to non-rectangular
shapes; and (b) it can be run in parallel using the `future` package.
This is particularly advantageous when using `winmove_agg`.
Parallelisation can be achieved by adding the following line prior to
running either `winmove_agg` or `nomove_agg`:

``` {tabindex="0"}
plan(multiprocess)
```

### Other features

In testing the method, we needed to simulate landscapes (which we did
with a great rOpenSci package
[`NLMR`](https://ropensci.github.io/NLMR/)). To avoid edge effects, we
wanted to create the effect of a torus (an infinite surface where the
top joins the bottom, and the left side meets the right side), but could
not find a function to do so - so we wrote one instead. Hopefully others
will find this useful when working with simulated landscapes.

::: {.highlight}
``` {.chroma tabindex="0"}
torus <- create_torus(cat_ls, 15)p1 <- show_landscape(cat_ls, legend_title = "", xlab = "", ylab = "") + ggtitle("Landscape")p2 <- show_landscape(torus, legend_title = "", xlab = "", ylab = "") + ggtitle("Landscape + torus")p1 + p2
```
:::

![](img/blog-images/2019-07-16-grainchanger/example2-1.png)

### Contributing

We hope to continue to build on the functions so that they are as
optimised as we can make them. If you have a particular function you
would like to see included in the package, please create a new
[issue](https://github.com/ropensci/grainchanger/issues), or submit a
pull request making sure develop is the destination branch on the
grainchanger repository.

In the long term, we would like to move the package into C++ so that all
functions are fast. If you are interested in collaborating on this
please get in touch.

### Acknowledgements

I would really like to acknowledge the [rOpenSci review
process](https://github.com/ropensci/software-review/issues/289). The
open nature of the reviews made for a thoughtful process, and I learned
a lot from it. Thanks to [Max Joseph](https://github.com/mbjoseph) and
[John Baumgartner](https://github.com/johnbaums) for their time and
expertise. The package and it's documentation are much improved thanks
to their input. [Felix
Eigenbrod](https://www.southampton.ac.uk/geography/about/staff/fe1y09.page)
and I worked on the methods behind the package together as part of the
ERC funded project *Scaling rules for ecosystem services*. Thanks also
for [Marco Sciaini](https://github.com/marcosci) for helpful discussion
around how to optimise some of this code and the tips on R package
development.

::: {.footnotes role="doc-endnotes"}

------------------------------------------------------------------------

1. ::: {#fn:1}
 [Graham et al. (2019) *Methods in Ecology and Evolution*.
 10:767--778](https://doi.org/10.1111/2041-210X.13177) [↩︎](#fnref:1){.footnote-backref
 role="doc-backlink"}
 :::
:::

The grainchanger package provides functionality for data aggregation to a coarser resolution via moving-window or direct methods. Why do we need new methods for data aggregation? As landscape ecologists and macroecologists, we often need to aggregate data in order to harmonise datasets. In doing so, we often lose a lot of information about the spatial structure and environmental heterogeneity of data measured at finer resolution.

Aggregating spatial data with the grainchanger package

Alec Robitaille

Quinn Webber

Eric Vander Wal

[spatsoc](https://github.com/ropensci/spatsoc) is an R package written
by Alec Robitaille, Quinn Webber and Eric Vander Wal of the Wildlife
Evolutionary Ecology Lab ([WEEL](https://weel.gitlab.io)) at Memorial
University of Newfoundland. It is the lab's first R package and was
recently accepted through the rOpenSci onboarding process with a big
thanks to reviewers [Priscilla Minotti](https://github.com/pmnatural)
and [Filipe Teixeira](https://github.com/FilipeamTeixeira), and editor
[Lincoln Mullen](https://github.com/lmullen).

`spatsoc` started as a single function (what would eventually become
`group_pts`) written by Alec in 2017 to help answer some of the
questions that Quinn and Eric were asking about how animal social
structure is related to spatial processes. These ideas were originally
proposed by Quinn and Eric in their recent review paper ^[1](#ref1)^.
After our ideas were published, we began using this early function to
determine when GPS collared caribou (*Rangifer tarandus*) in
Newfoundland were recorded within 50 m of one another, within 5 minutes.
This spatiotemporal grouping allowed us to build and analyze social
association networks with
[`asnipe`](https://sites.google.com/site/drfarine/r-packages/asnipe) and
[`igraph`](https://igraph.org/r/). Using animal telemetry data with
social network analysis allowed us to draw new insights from a long term
movement dataset.

<figure class="center">
<img
src="/img/blog-images/2018-12-04-spatsoc/Caribou_Juliana_Balluffi-Fry.jpg"
style="width:70.0%" alt=" Juliana Ballufi-Fry" />
<figcaption aria-hidden="true"><a
href="https://www.instagram.com/juliana.wildlife/"></a> <a
href="https://twitter.com/balluffi"></a>Juliana
Ballufi-Fry</figcaption>
</figure>

### Social network analysis

Social network analysis (SNA) is a powerful tool for measuring
associations and relationships between individuals with implications for
disease transmission, species interactions and community structure. SNA
uses graph theory and networks to reconstruct social structures, where
nodes represent individuals and edges represent the relationships or
associations between them. Social networks built from `spatsoc` results
are explicitly association networks, with the edges being defined by the
spatiotemporal grouping of relocations as opposed to observed social
interactions. Spiegel et al. ^[2](#ref2)^ showed that combining social
network analysis with telemetry data of sleepy lizards can reveal
insights about social structures in which the individuals and their
behaviours that shape it are otherwise difficult to directly observe.

### Animal telemetry data

The quantity and accessibility of animal movement data is increasing. A
number of recent articles have suggested we are in a "golden age of
animal tracking science" with large, extensive movement datasets being
collected for a diversity of species, many of which are publicly
available ^[3](#ref3)^. [Movebank](https://www.movebank.org/node/2) is a
global data archive of animal movement data, with over 800 taxa and 1
billion relocations. It also hosts an archive of published datasets with
accompanying DOI and licenses used in over 100 publications, backed by a
formal review process to ensure data quality.

Animal telemetry technology is improving rapidly, as battery systems
improve and sensors become smaller. Deployments of more advanced
biologging devices on animals means less frequent captures to redeploy
and more continuous data. [ICARUS](https://www.icarus.mpg.de/) (the
International Cooperation for Animal Research Using Space) is a
collaborative project using GPS and 3D accelerometry loggers powered by
solar panels and rechargeable batteries to produce continuous, high
temporal and spatial resolution data. The goal of ICARUS is to use this
data to better understand movement of animals and selecting forces that
act on individuals and contribute to our understanding of health,
disease epidemics and even potentially [predict volcanoes and
earthquakes](https://www.icarus.mpg.de/28810/animals-warning-sensors)!

> spatsoc allows us to answer the question: when and where are
> individuals close to each other?

The spatiotemporal grouping of `spatsoc` can lead to many outcomes. We
developed it to build proximity based social association networks from
GPS relocation data but others may use it to measure species
interactions or shared space use among con- or hetero-specifics. We are
not restricted to dyadic interactions since `spatsoc` groups relocations
for all individuals in space and time. spatsoc's flexible functions let
user define spatial and temporal grouping thresholds.

#### Package installation

Install the development version with
[`remotes`](https://github.com/r-lib/remotes).

::: {.highlight}
``` {.chroma tabindex="0"}
remotes::install_github('ropensci/spatsoc')# or CRANinstall.packages('spatsoc')
```
:::

#### Temporal and point based spatial grouping

For this example, we will use bison data from the "Dunn Ranch Bison
Tracking Project" (Max Planck Institute and Nature Conservancy Missouri)
available on Movebank
[here](https://www.movebank.org/panel_embedded_movebank_webapp?gwt_fragment=page%3Dsearch_map_linked%2CindividualIds%3D8070814%2B8070814%2B8070814%2B8070774%2B8070774%2B8070774%2B8070775%2B8070813%2B8070813%2B8070813%2B8070749%2B8057017%2B8070744%2B8070744%2B8070744%2B54137325%2B54137415%2B54137359%2B13848400%2B13848400%2B13848400%2B13848400%2B8070773%2B8070812%2B8070812%2B8070735%2B8070811%2Clat%3D40.45725700561854%2Clon%3D-94.10498513423153%2Cz%3D13).
We will import the downloaded data (CSV format) with
[`data.table`](https://github.com/Rdatatable/data.table)'s `fread`. This
function is fast when data gets big and `data.table` is the required
format for all `DT` arguments in spatsoc.

Note: to download the data from Movebank you will need to make an
account and accept the terms of the data. We are only going to use a
subset of the available data for this example and we have renamed the
downloaded file to `Bison.csv` for brevity. Lastly, select "Add UTM
coordinates" on the download page to ensure we have projected
coordinates ready for spatsoc.

::: {.highlight}
``` {.chroma tabindex="0"}
# Load packageslibrary(spatsoc)library(data.table)# Read data as a data.tableDT <- fread('Bison.csv')# Cast the character column 'timestamp' as a POSIXct, saved as 'datetime'# GMT/UTC timezone as specific in the data's README.txtDT[, datetime := as.POSIXct(timestamp, tz = 'GMT')]# Subset to 2015DT <- DT[year(datetime) == 2015]# Temporal grouping with threshold of 5 minutesgroup_times(DT = DT, datetime = 'datetime', threshold = '5 minutes')# Spatial grouping with threshold of 25 metersgroup_pts(DT = DT, threshold = 25, id = 'individual_local_identifier', coords = c('utm-easting', 'utm-northing'), timegroup = 'timegroup')
```
:::

Behind the scenes, we are building a distance matrix and using the
connected components to determine group membership.

::: {.highlight}
``` {.chroma tabindex="0"}
library(igraph)distMatrix <- as.matrix(dist(cbind(X, Y), method = "euclidean"))graphAdj <- graph_from_adjacency_matrix(distMatrix <= threshold, mode = 'undirected', diag = FALSE)groups <- clusters(graphAdj)$membership
```
:::

A distance matrix measures the spatial distance between all of the
individuals. The spatial distance `threshold` defined by the user is
used to turn the distance matrix into a boolean matrix.

<figure class="center">
<img src="/img/blog-images/2018-12-04-spatsoc/dist-matrix.png" />
</figure>

This matrix is then passed to `igraph` to build a graph, and the
connected components are extracted. The connected components represent
the distinct, separate subgraphs that make up a graph. In this case,
there are three connected components representing our spatiotemporal
groups.

<figure class="center">
<img
src="/img/blog-images/2018-12-04-spatsoc/connected-components.png" />
</figure>

### Reflecting on the process

`spatsoc` was submitted to rOpenSci on July 16 2018 and accepted on
October 26 2018. Based on feedback during the review process, the
package documentation and usability were greatly improved. The rOpenSci
onboarding process represented a different standard than CRAN, where
test coverage, useful documentation and
[`goodpractice`](https://github.com/MangoTheCat/goodpractice) were just
as important as getting `0 errors | 0 warnings | 0 messages`. We really
appreciated the common goal of improving the package shared by package
authors, rOpenSci reviewers, and editor.

The package building process gave us a better understanding of how the
functions and packages we use in our analyses are built. Using powerful
tools like unit testing and CRAN checks we could ensure our code was
well written. Alec was really surprised how much he enjoyed writing unit
tests, thinking of all the different ways things could fail and setting
up the necessary framework to prevent it, returning useful warnings and
errors when it did.

In the future: more packages! For personal use, for fun, and for sharing
new methods and tools with the R community. The flexibility and power of
the package as a fundamental unit of sharing functions (and data,
projects, etc.). Next time, we will write unit tests and documentation
as we go (hopefully!). We are looking forward to engaging in more
collaboration opportunities in the R community through the rOpenSci
onboarding process as authors and reviewers or less formally through
Pull Requests or Issue Boards.

### Resources

Check out the [repository](https://github.com/ropensci/spatsoc),
contribute to [open issues](https://github.com/ropensci/spatsoc/issues)
or [open a new issue](https://github.com/ropensci/spatsoc/issues/new)!
We have written three vignettes and created a
[website](https://spatsoc.gitlab.io) using
[pkgdown](https://pkgdown.r-lib.org/). If the documentation and these
resources still leave you with questions about using the package, please
contact us or create an issue.

We are always open to chat about our ideas, the research we are
currently doing with spatsoc, and the research we plan to do in the
future, so please let us know if you have any questions.

### References

[1]{#ref1}: Webber, Q. M., & Vander Wal, E. (2018). [An evolutionary
framework outlining the integration of individual social and spatial
ecology.](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/1365-2656.12773)
*Journal of Animal Ecology*.

[2]{#ref2}: Spiegel, O., Sih, A., Leu, S. T., & Bull, C. M. (2018).
[Where should we meet? Mapping social network interactions of sleepy
lizards shows sex-dependent social network
structure.](https://www.sciencedirect.com/science/article/pii/S0003347217303639)
*Animal Behaviour*.

[3]{#ref3}: Kays, R., Crofoot, M. C., Jetz, W., & Wikelski, M. (2015).
[Terrestrial animal tracking as an eye on life and
planet.](https://science.sciencemag.org/content/348/6240/aaa2478)
*Science*.

spatsoc is an R package written by Alec Robitaille, Quinn Webber and Eric Vander Wal of the Wildlife Evolutionary Ecology Lab (WEEL) at Memorial University of Newfoundland. It is the lab’s first R package and was recently accepted through the rOpenSci onboarding process with a big thanks to reviewers Priscilla Minotti and Filipe Teixeira, and editor Lincoln Mullen.

Messaggi di Rogue Scholar

Working with Qualtrics Data - Part 2: Excluding Data

coder Makes Medical Coding less Messy

Introducing BaseSet for mathematical sets

Aggregating spatial data with the grainchanger package

Detecting spatiotemporal groups in relocation data with spatsoc