Added several things

Most notably documentation of remaining functions located in helpfulfunctions.R

.Rbuildignore

@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,45 @@
+# History files
+.Rhistory
+.Rapp.history
+
+# Session Data files
+.RData
+
+# User-specific files
+.Ruserdata
+
+# Example code in package build process
+*-Ex.R
+
+# Output files from R CMD build
+/*.tar.gz
+
+# Output files from R CMD check
+/*.Rcheck/
+
+# RStudio files
+.Rproj.user/
+
+# produced vignettes
+vignettes/*.html
+vignettes/*.pdf
+
+# OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
+.httr-oauth
+
+# knitr and R markdown default cache directories
+*_cache/
+/cache/
+
+# Temporary files created by R markdown
+*.utf8.md
+*.knit.md
+
+# R Environment Variables
+.Renviron
+
+# pkgdown site
+docs/
+
+# translation temp files
+po/*~

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: bremla
 Type: Package
 Title: Bayesian Regression Modeling of Layered Archives
-Version: 0.1.4
+Version: 0.1.5
 Author: Eirik Myrvoll-Nilsen
 Maintainer: Eirik Myrvoll-Nilsen <[email protected]>
 Description: Performs efficient Bayesian regression modeling of layer-counted climate records. Default method of inference is integrated nested Laplace approximations (INLAs)

NAMESPACE

@@ -7,15 +7,20 @@ export(bremla_biased_chronologies)
 export(bremla_chronology_simulation)
 export(bremla_modelfitter)
 export(bremla_prepare)
+export(bremla_simulationsummarizer)
+export(linramp)
 export(linrampfitter)
+export(meanmaker)
 export(rgeneric.uneven.AR1)
+export(which.index)
 import(INLA)
 import(Matrix)
 import(matrixStats)
 import(numDeriv)
 importFrom(INLA,inla)
 importFrom(INLA,inla.ar.pacf2phi)
 importFrom(INLA,inla.emarginal)
+importFrom(INLA,inla.hpdmarginal)
 importFrom(INLA,inla.hyperpar.sample)
 importFrom(INLA,inla.rgeneric.define)
 importFrom(INLA,inla.smarginal)

R/bremla_prepare.R

@@ -47,24 +47,30 @@ bremla_prepare = function(age,depth,proxy, events=NULL,nsims=10000, eventmeasure
  if(!is.null(events)){
  eventindexes = numeric(length(events))
 
- for(i in 1:length(events)){
- if(tolower(eventmeasure) %in% c("depth","z")){
- if(events[i] < min(z) || events[i]>max(z)){
- warning(paste0("Event ",i,", located at ",events[i]," is outside the interval covered by 'depth' (",min(z),", ",max(z),"). The event will be omitted!"))
- eventindexes[i] = NA
- }else{
- eventindexes[i] = which(abs(events[i]-z) == min(abs(events[i]-z)))
- }
- }else if(tolower(eventmeasure) %in% c("age","time","y")){
- if(events[i] < min(y) || events[i]>max(y)){
- warning(paste0("Event ",i,", located at ",events[i]," is outside the interval covered by 'age' (",min(age),", ",max(age),"). The event will be omitted!"))
- eventindexes[i] = NA
- }else{
- eventindexes[i] = which(abs(events[i]-y) == min(abs(events[i]-y)))
- }
- }
-
+ if(tolower(eventmeasure) %in% c("depth","z")){
+ eventindexes = which.index (events, z)
+ }else if(tolower(eventmeasure) %in% c("age","time","y")){
+ eventindexes = which.index (events, y)
  }
+ # for(i in 1:length(events)){
+ # if(tolower(eventmeasure) %in% c("depth","z")){
+ # eventindexes = which.index (events, record)
+ # if(events[i] < min(z) || events[i]>max(z)){
+ # warning(paste0("Event ",i,", located at ",events[i]," is outside the interval covered by 'depth' (",min(z),", ",max(z),"). The event will be omitted!"))
+ # eventindexes[i] = NA
+ # }else{
+ # eventindexes[i] = which(abs(events[i]-z) == min(abs(events[i]-z)))
+ # }
+ # }else if(tolower(eventmeasure) %in% c("age","time","y")){
+ # if(events[i] < min(y) || events[i]>max(y)){
+ # warning(paste0("Event ",i,", located at ",events[i]," is outside the interval covered by 'age' (",min(age),", ",max(age),"). The event will be omitted!"))
+ # eventindexes[i] = NA
+ # }else{
+ # eventindexes[i] = which(abs(events[i]-y) == min(abs(events[i]-y)))
+ # }
+ # }
+ #
+ # }
  eventindexes = unique(c(1,eventindexes[!is.na(eventindexes)]))
  nevents = length(eventindexes)
  }

R/helpfulfunctions.R

@@ -1,3 +1,18 @@
+#' Mean vector from fixed effects
+#'
+#' Gives the mean number of layer differences per depth from the fixed effects.
+#'
+#' @param coefs Fixed effects
+#' @param reg.model list specifying the structure of the linear regression model.
+#' @param nevents The number of climate transitions
+#' @param data data.frame obtained from \code{\link{bremla_prepare}}.
+#'
+#' @return Returns a vector representing the mean number of layer difference per depth interval.
+#' @author Eirik Myrvoll-Nilsen, \email{[email protected]}
+#' @seealso \code{\link{bremla_prepare},\link{bremla_chronology_simulation}}
+#' @keywords bremla mean
+#'
+#' @examples
 #' @export
 meanmaker = function(coefs,reg.model,nevents=69,data){
  ## Computes mean vector from given fixed effects 'coefs'.
@@ -36,6 +51,22 @@ meanmaker = function(coefs,reg.model,nevents=69,data){
  }
  return(fitted)
 }
+#' Linear ramp function
+#'
+#' Computes the (noiseless) linear ramp function.
+#'
+#' @param t Vector describing the x-axis (here, depth)
+#' @param t0 the onset of the transition
+#' @param dt The transition duration
+#' @param y0 The initial level of the observed values (here, proxy values)
+#' @param dy The increase in the observed values
+#'
+#' @return Returns a vector y-axis corresponding to \code{t}
+#' @author Eirik Myrvoll-Nilsen, \email{[email protected]}
+#' @seealso \code{\link{linrampfitter}}
+#' @keywords linramp
+#'
+#' @examples
 #' @export
 linramp = function(t,t0=0,dt=1,y0=0,dy=1){
  y = numeric(length(t))
@@ -57,6 +88,19 @@ linramprev = function(t,t0=0,dt=1,y0=0,dy=1){
 }
 
 
+#' Corresponding index
+#'
+#' Finds the indices where events best match values in a given vector.
+#'
+#' @param events Vector describing the values of interest (here: either depth or age of events)
+#' @param record Vector we want to search (here: either full depth or age record)
+#'
+#' @return Returns a vector of indices for which elements of \code{record} best matches the elements of \code{events}-
+#' @author Eirik Myrvoll-Nilsen, \email{[email protected]}
+#' @seealso \code{\link{bremla}}
+#' @keywords indices
+#'
+#' @examples
 #' @export
 which.index = function(events, record){ ## Finds which indices of 'record' that are located closest to the 'event's
  eventindexes = numeric(length(events))
@@ -69,13 +113,28 @@ which.index = function(events, record){ ## Finds which indices of 'record' that
  }
 
  }
- eventindexes = unique(c(1,eventindexes[!is.na(eventindexes)])) #Placing transition at the start of record. Removing NA and duplicates
+ #eventindexes = unique(c(1,eventindexes[!is.na(eventindexes)])) #Placing transition at the start of record. Removing NA and duplicates
  return(eventindexes)
 }
+
+#' Simulation-summarizer
+#'
+#' Computes posterior marginal mean and uncertainty intervals from simulations.
+#'
+#' @param object List object which is the output of function \code{\link{bremla_chronology_simulation}}
+#' @param interval character describing which uncertainty intervals should be used. Highest posterior density is used by default
+#' @param print.progress Boolean describing whether or not progress should be printed on screen.
+#'
+#' @return Returns the \code{object} list from the input and appends additional summary statistics: posterior marginal mean and uncertainty interval-
+#' @author Eirik Myrvoll-Nilsen, \email{[email protected]}
+#' @seealso \code{\link{bremla_chronology_simulation}}
+#' @keywords indices
+#'
+#' @examples
 #' @export
 #' @importFrom matrixStats rowMeans2 rowVars
 #' @importFrom INLA inla.hpdmarginal
-bremla_simulationsummarizer = function(object,print.progress=FALSE){
+bremla_simulationsummarizer = function(object,interval="hpd",print.progress=FALSE){
  if(print.progress) cat("Computing posterior marginal mean and 95% hpd intervals from chronology samples...\n",sep="")
  time.start = Sys.time()
  n = dim(object$simulation$age)[1]
@@ -89,8 +148,8 @@ bremla_simulationsummarizer = function(object,print.progress=FALSE){
  hpdupper[i] = inla.hpdmarginal(0.95,dens)[2]
  }
  time.summary = Sys.time()
+ object$simulation$summary = list(mean=meanvek,sd=sdvek,hpd0.025=hpdlower,hpd0.975=hpdupper, sim.sum.time=time.summary,.args=interval=interval,print.progress=print.progress)
  if(print.progress) cat(" completed in ",difftime(time.summary,time.start,units="secs")[[1]],"\n",sep="")
- object$simulation$summary = list(mean=meanvek,sd=sdvek,hpd0.025=hpdlower,hpd0.975=hpdupper, sim.sum.time=time.summary)
  object$time$samplesummary = list(total=difftime(time.summary,time.start,units="secs")[[1]])
  return(object)
 }

bremla.Rproj

@@ -0,0 +1,20 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source