Skip to content

Commit

Permalink
fixed initial value problem
Browse files Browse the repository at this point in the history
  • Loading branch information
@abigailkeller
abigailkeller committed Sep 13, 2024
1 parent e3d9752 commit fc577af
Show file tree
Hide file tree
Showing 3 changed files with 271 additions and 146 deletions.
158 changes: 79 additions & 79 deletions .Rhistory
View file
Original file line number Diff line number Diff line change
@@ -1,82 +1,3 @@
paste("Initial values for 'mu' should be numeric values > 0.",
paste0('See the eDNAjoint guide for help formatting ',
'initial values: '),
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase1.html#initialvalues'),
sep='\n'))
#38. initial values check: mu length
n.chain <- 4
inits <- list()
for(i in 1:n.chain){
inits[[i]] <- list(
mu <- stats::runif(4, 0.1, 1),
p10 <- stats::runif(1,log(0.0001),log(0.08)),
alpha <- rep(0.1,3)
)
names(inits[[i]]) <- c('mu','p10','alpha')
}
site.cov <- rbind(c(4,1),c(1,1))
colnames(site.cov) <- c('var_a','var_b')
expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1),c(1,1,NA)),
qPCR.N = rbind(c(3,3,3),c(3,3,NA)),
count = rbind(c(4,1,1),c(1,1,NA)),
site.cov = site.cov),
cov = c('var_a','var_b'), initial_values = inits,
multicore = FALSE),
paste(paste0("The length of initial values for 'mu' should ",
"equal the number of sites."),
paste0('See the eDNAjoint guide for help formatting ',
'initial values: '),
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase1.html#initialvalues'),
sep='\n'))
#39. initial values check: p10 numeric
n.chain <- 4
inits <- list()
for(i in 1:n.chain){
inits[[i]] <- list(
mu <- stats::runif(2,0,1),
p10 <- "-1",
alpha <- rep(0.1,3)
)
names(inits[[i]]) <- c('mu','p10','alpha')
}
site.cov <- rbind(c(4,1),c(1,1))
colnames(site.cov) <- c('var_a','var_b')
expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1),c(1,1,NA)),
qPCR.N = rbind(c(3,3,3),c(3,3,NA)),
count = rbind(c(4,1,1),c(1,1,NA)),
site.cov = site.cov),
cov = c('var_a','var_b'), initial_values = inits,
multicore = FALSE),
paste("Initial values for 'p10' should be numeric.",
paste0('See the eDNAjoint guide for help formatting ',
'initial values: '),
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase1.html#initialvalues'),
sep='\n'))
#40. initial values check: p10 length
n.chain <- 4
inits <- list()
for(i in 1:n.chain){
inits[[i]] <- list(
mu <- stats::runif(2,0,1),
p10 <- stats::runif(2,log(0.0001),log(0.08)),
alpha <- rep(0.1,3)
)
names(inits[[i]]) <- c('mu','p10','alpha')
}
site.cov <- rbind(c(4,1),c(1,1))
colnames(site.cov) <- c('var_a','var_b')
expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1),c(1,1,NA)),
qPCR.N = rbind(c(3,3,3),c(3,3,NA)),
count = rbind(c(4,1,1),c(1,1,NA)),
site.cov = site.cov),
cov = c('var_a','var_b'), initial_values = inits,
multicore = FALSE),
paste("The length of initial values for 'p10' should equal 1.",
paste0('See the eDNAjoint guide for help formatting ',
'initial values: '),
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase1.html#initialvalues'),
sep='\n'))
Expand Down Expand Up @@ -510,3 +431,82 @@ color = 'survey type')+
ggplot2::scale_color_manual(values = scales::hue_pal()(length(names)),
labels = names)+
ggplot2::theme_minimal()
roxygen2::roxygenise()
install.packages('srr')
options(repos = c(
ropenscireviewtools = "https://ropensci-review-tools.r-universe.dev",
CRAN = "https://cloud.r-project.org"))
install.packages("srr")
roxygen2::roxygenise()
GLNA_grab <- read.csv('../neha_data/GLNA_grab.csv')
OPAC_grab <- read.csv('../neha_data/OPAC_grab.csv')
PRLI_grab <- read.csv('../neha_data/PRLI_grab.csv')
GLNA_dna <- read.csv('../neha_data/eDNAJoint_GLNA_B.csv')
OPCA_dna <- read.csv('../neha_data/eDNAJoint_OPAC_BD.csv')
PRLI_dna <- read.csv('../neha_data/eDNAJoint_PRLI_AW.csv')
OPCA_grab <- read.csv('../neha_data/OPAC_grab.csv')
GLNA_grab <- read.csv('../neha_data/GLNA_grab.csv')
OPCA_grab <- read.csv('../neha_data/OPAC_grab.csv')
PRLI_grab <- read.csv('../neha_data/PRLI_grab.csv')
GLNA_dna <- read.csv('../neha_data/eDNAJoint_GLNA_B.csv')
OPCA_dna <- read.csv('../neha_data/eDNAJoint_OPAC_BD.csv')
PRLI_dna <- read.csv('../neha_data/eDNAJoint_PRLI_AW.csv')
qPCR_N_wide <- OPCA_dna %>%
pivot_wider(id_cols=site,
names_from=biological_replicate,
values_from=N) %>%
select(-site)
# K (number of positive qPCR detections per biological replicate (i.e., water sample))
qPCR_K_wide <- OPCA_dna %>%
pivot_wider(id_cols=site,
names_from=biological_replicate,
values_from=K) %>%
select(-site)
library(tidyverse)
qPCR_N_wide <- OPCA_dna %>%
pivot_wider(id_cols=site,
names_from=biological_replicate,
values_from=N) %>%
select(-site)
# K (number of positive qPCR detections per biological replicate (i.e., water sample))
qPCR_K_wide <- OPCA_dna %>%
pivot_wider(id_cols=site,
names_from=biological_replicate,
values_from=K) %>%
select(-site)
# counts (number of individuals in traditional survey sample)
count_wide <- OPCA_grab %>%
pivot_wider(id_cols=site,
names_from=grab_replicate,
values_from=count) %>%
select(-site)
plotdata <- data.frame(
x = rowMeans(count_wide,na.rm=TRUE),
y = rowMeans(qPCR_K_wide,na.rm=TRUE),
label = 1:dim(count_wide)[1]
)
ggplot(data=plotdata,aes(x=x,y=y))+
geom_point()+
geom_text(aes(label=factor(label), vjust = -1))+
#scale_y_continuous(limits=c(3.5,6.8))+
labs(x='mean grab count',y='mean qPCR detections')+
ggtitle('GLNA')
ggplot(data=plotdata,aes(x=x,y=y))+
geom_point()+
geom_text(aes(label=factor(label), vjust = -1))+
#scale_y_continuous(limits=c(3.5,6.8))+
labs(x='mean grab count',y='mean qPCR detections')+
ggtitle('OPCA')
# combine all data into one named list
# note: make sure to convert the wide dataframes into matrices
data <- list(
qPCR.N = as.matrix(qPCR_N_wide),
qPCR.K = as.matrix(qPCR_K_wide),
count = as.matrix(count_wide)
)
# OPCA grab
fit <- jointModel(data=data, family='negbin', p10priors = c(1,10))
library(eDNAjoint)
# OPCA grab
fit <- jointModel(data=data, family='negbin', p10priors = c(1,10))
jointSummarize(fit$model)
123 changes: 123 additions & 0 deletions .ipynb_checkpoints/README-checkpoint.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,123 @@

<!-- README.md is generated from README.Rmd. Please edit that file -->

# eDNAjoint

<img src="man/figures/logo.png" align="right" height="200" dpi="700"/>

<!-- badges: start -->

[![R-CMD-check](https://github.com/abigailkeller/eDNAjoint/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/abigailkeller/eDNAjoint/actions/workflows/R-CMD-check.yaml)
[![codecov](https://codecov.io/gh/abigailkeller/eDNAjoint/graph/badge.svg?token=AEVR9NSQ9Z)](https://codecov.io/gh/abigailkeller/eDNAjoint)
[![Status at rOpenSci Software Peer
Review](https://badges.ropensci.org/642_status.svg)](https://github.com/ropensci/software-review/issues/642)
[![Lifecycle:
stable](https://img.shields.io/badge/lifecycle-stable-green.svg)](https://lifecycle.r-lib.org/articles/stages.html#stable)

<!-- badges: end -->

The package *eDNAjoint* is useful for interpreting observations from
paired or semi-paired environmental DNA (eDNA) and traditional surveys.
The package runs a Bayesian model that integrates these two data streams
to jointly estimate parameters like the false positive probability of
eDNA detection and expected catch rate at a site. Optional model
variations allow inclusion of site-level covariates that scale the
sensitivity of eDNA sampling relative to traditional sampling, as well
as estimation of catchability coefficients when multiple traditional
gear types are used. Additional functions in the package facilitate
interpretation of model fits.

<div style="text-align: center;">

<img src="man/figures/basic_diagram_final-01.png" height="400" dpi="700"/>

</div>

## Installation

You can install the development version of eDNAjoint. Note that
installation time from source can be high, as unlike binary file formats
provided by Cran, the program needs to compile during installation:

``` r
library(devtools)
install_github('abigailkeller/eDNAjoint')
```

## Example

The main functionality in *eDNAjoint* is the use of `jointModel()` that
will fit the model to data. Further functions like `jointSummarize()`
and `detectionCalculate()` can be used to help with model fit
interpretation.

This example fits the joint model to data from paired, replicated eDNA
qPCR and seine sampling observations of endangered tidewater gobies
(*Eucyclogobius newberryi*) from a study by Schmelzle and Kinziger
(2016). The following variation of the joint model includes site-level
covariates that scale the sensitivity of eDNA sampling relative to
traditional sampling.

``` r
library(eDNAjoint)
data(gobyData)
# run the joint model with two covariates
goby.fit <- jointModel(data = gobyData, cov = c('Filter_time','Salinity'),
family = 'poisson', p10priors = c(1,20), q = FALSE)
```

And then this model fit can be accessed to do things like summarize the
posterior distribution for the probability of a false positive
detection, $p_{10}$:

``` r
# summarize p10 posterior
jointSummarize(goby.fit$model, par = 'p10')
#> mean se_mean sd 2.5% 97.5% n_eff Rhat
#> p10 0.003 0 0.001 0.001 0.007 15361.8 1
```

Or to find the number of eDNA samples and traditional survey samples
necessary to detect presence of the species at a given expected catch
rate:

``` r
# find the number of samples necessary to detect presence with 0.9 probability at the mean covariate values,
# if the expected catch rate (mu) is 0.1, 0.5, or 1 individuals/traditional survey unit.
detectionCalculate(goby.fit$model, mu = c(0.1,0.5,1),
cov.val = c(0,0), probability = 0.9)
#> mu n_traditional n_eDNA
#> [1,] 0.1 24 14
#> [2,] 0.5 5 4
#> [3,] 1.0 3 2
```

## Vignette

You can find much more detailed examples of the functions in *eDNAjoint*
and the model underlying the package in the [package
vignette](https://bookdown.org/abigailkeller/eDNAjoint_vignette/).

## Citation

Keller, A.G. (2024). *eDNAjoint: R package for interpreting paired and
semi-paired environmental DNA and traditional surveys*.
<https://github.com/abigailkeller/eDNAjoint>.

## Code of Conduct

Please note that eDNAjoint is released with a [Contributor Code of
Conduct](https://ropensci.org/code-of-conduct/). By contributing to this
project you agree to abide by its terms.

## References

Keller, A.G., Grason, E.W., McDonald, P.S., Ramon-Laca, A., Kelly, R.P.
(2022). Tracking an invasion front with environmental DNA. *Ecological
Applications*. 32(4): e2561. <https://doi.org/10.1002/eap.2561>

Schmelzle, M.C. and Kinziger, A.P. (2016). Using occupancy modelling to
compare environmental DNA to traditional field methods for
regional-scale monitoring of an endangered aquatic species. *Molecular
Ecology Resources*. 16(4): 895-908.
<https://doi.org/10.1111/1755-0998.12501>
Loading

0 comments on commit fc577af

Please sign in to comment.