Local Polynomial Regression of State Dependent Covariates in State-Space Models

Facilitates local polynomial regression for state dependent covariates in state-space models. The functionality can also be used from 'C++' based model builder tools such as 'Rcpp'/'inline', 'TMB', or 'JAGS'.


library(knitr)
opts_knit$set(root.dir=normalizePath('./'))
opts_chunk$set(fig.path = "./tools/README_fig/", dev='png')

Installing

covafillr version can be installed from CRAN with

install.packages("covafillr")

The development version of covafillr can be installed with

devtools::install_github("calbertsen/covafillr")

Installing with JAGS

To use covafillr with JAGS, JAGS must be installed before the covafillr package, and the package must be installed using the same compiler as JAGS is installed with.

On Unix(-like) systems, pkg-config is used to find the relevant paths to compile covafillr against JAGS, such as

pkg-config --cflags jags
pkg-config --libs jags
## -I/usr/local/include/JAGS 
## -L/usr/local/lib -ljags

The package will only be installed with the JAGS module if the configure argument --with-jags is given. Note that the package must be compiled with the same compiler as JAGS.

On Windows, the R package rjags is used to find the paths. covafillr can be installed without using rjags by setting a system variable JAGS_ROOT to the folder where JAGS is installed, e.g., by running

Sys.setenv(JAGS_ROOT='C:/Program Files/JAGS/JAGS-4.1.0')

before installation. Similar to Unix(-like) systems, the package is only installed with the JAGS module if the system variable USE_JAGS is set, e.g., by running

Sys.setenv(USE_JAGS='yes')

Simple examples

Note that more examples are available in the inst/examples folder.

Using from R

library(methods)
library(covafillr)

The package can be used from R to do local polynomial regression and a search tree approximation of local polynomial regression. Both are implemented with reference classes.

Local polynomial regression

The reference class for local polynomial regression is called covafill.

methods::getRefClass('covafill')
## Generator for class "covafill":
## 
## Class fields:
##                   
## Name:          ptr
## Class: externalptr
## Locked Fields: "ptr"
## 
## 
## Class Methods: 
##      "import", "getDegree", ".objectParent", "setBandwith", "residuals", 
##      "usingMethods", "show", "getClass", "untrace", "export", 
##      "copy#envRefClass", "initialize", ".objectPackage", "callSuper", 
##      "getDim", "copy", "getBandwith", "predict", "initFields", 
##      "getRefClass", "trace", "field"
## 
## Reference Superclasses: 
##      "envRefClass"

To illustrate the usage we simulate data.

fn <- function(x) x ^ 4 - x ^ 2
x <- runif(2000,-10,10)
y <- fn(x) + rnorm(2000,0,0.1)

An object of the reference class is created by

cf <- covafill(coord = x,obs = y,p = 5L)

where p is the polynomial degree. The bandwith can be set by the argument h. By default, a value is suggested by the function suggestBandwith. Information about the class can be extracted (and changed) by the following functions:

cf$getDim()
## [1] 1
cf$getDegree()
## [1] 5
cf$getBandwith()
## [1] 15.2975
cf$setBandwith(1.0)
## [1] 1
cf$getBandwith()
## [1] 1

To do local polynomial regression at a point, the $predict function is used.

x0 <- seq(-3,3,0.1)
y0 <- cf$predict(x0)

The function returns a matrix of estimated function values and derivatives.

par(mfrow=c(3,1))
plot(x0,y0[,1], main = "Function")
lines(x0,fn(x0))
plot(x0, y0[,2], main = "First derivative")
lines(x0, 4 * x0 ^ 3 - 2 * x0)
plot(x0, y0[,3], main = "Second derivative")
lines(x0, 3 * 4 * x0 ^ 2 - 2)

Search tree approximation

The reference class for search tree approximation to local polynomial regression is covatree

methods::getRefClass('covatree')
## Generator for class "covatree":
## 
## Class fields:
##                   
## Name:          ptr
## Class: externalptr
## Locked Fields: "ptr"
## 
## 
## Class Methods: 
##      "import", ".objectParent", "usingMethods", "show", "getClass", 
##      "untrace", "export", "copy#envRefClass", "initialize", 
##      ".objectPackage", "callSuper", "getDim", "copy", "predict", 
##      "initFields", "getRefClass", "trace", "field"
## 
## Reference Superclasses: 
##      "envRefClass"

covatree has an aditional argument, minLeft, which is the minimum number of observations at which a sub tree will be created. Otherwise the functionality is similar.

ct <- covatree(coord = x,obs = y,p = 5L, minLeft = 50)
ct$getDim()
## [1] 1
y1 <- ct$predict(x0)
par(mfrow=c(2,1))
plot(x0,y1[,1], main = "Function")
lines(x0,fn(x0))
plot(x0, y1[,2], main = "First derivative")
lines(x0, 4 * x0 ^ 3 - 2 * x0)

Using with Rcpp/inline

The covafillr package provides a plugin for inline.

library(inline)

The following code does local polynomial regression at the point x based on the observations obs at the points coord. For convenience, the plugin provides the type definitions cVector and cMatrix to pass to the covafill constructor and operator.

cftest <- '
  cVector x1 = as<cVector>(x);
  cMatrix coord1 = as<cMatrix>(coord);
  cVector obs1 = as<cVector>(obs);
  int p1 = as<int>(p);
  cVector h1 = as<cVector>(h);
  covafill<double> cf(coord1,obs1,h1,p1);
  return wrap(cf(x1));'

An R function can now be defined with inlined C++ code using plugin='covafillr'.

fun <- cxxfunction(signature(x='numeric',
                             coord = 'matrix',
                             obs = 'numeric',
                             p = 'integer',
                             h = 'numeric'),
                   body = cftest,
                   plugin = 'covafillr'
                   )    
fun(c(0),matrix(x,ncol=1),y,2,1.0)
## [1] -0.03907458 -0.01879995

Using with TMB

To use covafillr with TMB, include covafill/TMB in the beginning of the cpp file.

// tmb_covafill.cpp  
#include <TMB.hpp>
#include <covafill/TMB>
 
template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_MATRIX(coord);
  DATA_VECTOR(covObs);
  DATA_INTEGER(p);
  DATA_VECTOR(h);
  PARAMETER_VECTOR(x);
  covafill<Type> cf(coord,covObs,h,p);
  Type val = evalFill((CppAD::vector<Type>)x, cf)[0];
  return val;
}

Instead of calling the usual operator, the covafill object is evaluated at a point with the evalFill function. This function enables TMB to use the estimated gradient in the automatic differentiation.

From R, the cpp file must be compiled with additional flags as seen below.

library(TMB)
TMB::compile('tmb_covafill.cpp',CXXFLAGS = cxxFlags())
## Note: Using Makevars in /home/cmoe/.R/Makevars

## [1] 0
dyn.load(dynlib("tmb_covafill"))

Then TMB can be used as usual.

dat <- list(coord = matrix(x,ncol=1),
            covObs = y,
            p = 2,
            h = 1.0)
 
obj <- MakeADFun(data = dat,
                 parameters = list(x = c(0)),
                 DLL = "tmb_covafill")
obj$fn(c(3.2))
## [1] 94.50085
obj$fn(c(0))
## [1] -0.03907458
obj$fn(c(-1))
## [1] -0.05840993
obj$gr()
## outer mgc:  3.661488

##           [,1]
## [1,] -3.661488

Using with rjags

library(rjags)
## Loading required package: coda

## Linked to JAGS 4.0.0

## Loaded modules: basemod,bugs

If covafillr is installed with JAGS, a module is compiled with the package. The module can be loaded in R by the function loadJAGSModule, a wrapper for rjags::load.module.

loadJAGSModule()
## module covafillr loaded

Then the function covafill is available to use in the JAGS code

# covafill.jags
model {
  cf <- covafill(x,obsC,obs,h,2.0)
  sigma ~ dunif(0,100)
  tau <- pow(sigma, -2)
  for(i in 1:N) {
    y[i] ~ dnorm(cf[i],tau)
  }
}
fun <- function(x) x ^ 2
n <- 100
x <- runif(n,-2,2)
y <- rnorm(n,fun(x),0.5)
obsC <- seq(-3,3,len=1000)
obs <- fun(obsC) + rnorm(length(obsC),0,0.1)

Then rjags can be used as usual.

jags <- jags.model('covafill.jags',
                   data = list(N = n,
                               x = matrix(x,ncol=1),
                               y = y,
                               obsC = matrix(obsC,ncol=1),
                               obs = obs,
                               h = c(1)),
                   n.chains = 1,
                   n.adapt = 100)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 1
##    Total graph size: 2314
## 
## Initializing model
samp <- jags.samples(jags,c('sigma','cf'),n.iter = 1000, thin = 10)
par(mfrow=c(2,1))
plot(x,samp$cf[,1,1])
hist(samp$sigma)

News

Covafillr News

CHANGES IN VERSION 0.4.2 BUG FIXES: - Rcpp:::Rcpp.plugin.maker changed to Rcpp::Rcpp.plugin.maker - Fixed index bug in tricubic interpolation

CHANGES IN VERSION 0.4.1

NEW FEATURES:
    - covafill and covatree now have $copy methods for non-shallow copying

BUG FIXES:
    - Files needed for README.md is now included in ./tools

OTHER CHANGES:
    - $show is now more informative for covafill and covatree reference classes

CHANGES IN VERSION 0.4.0

NEW FEATURES:
    - added stat_covafill to be used with the ggplot2 package
    - added kde function for kernel density estimates
    - added suggestBandwith function

OTHER CHANGES:
    - Changed calculations of standard deviations to follow Fan and Gijbels (1995)
    - Now possible to have polynomial degree < 1
    - Changed default polynomial degree to 3
    - Small changes to C++: Added destructors to tree methods; now registers routines; fixed pedantic warnings

CHANGES IN VERSION 0.3.0

NEW FEATURES:
    - It is now possible to get standard deviations for estimates with covafill.

BUG FIXES:
    - covafill$predict now set colnames correctly for p=2.
    - Weight matrix in covafill used sqrt of kernel weights.

OTHER CHANGES:
    - Installation with JAGS must now be explicitly invoked
    - $copy has been disabled for covafill and covatree reference classes.
    - .Object has been removed from argument list in $initialize  for covafill and covatree reference classes.

VERSION 0.2.1: First publicly available version

Reference manual

It appears you don't have a PDF plugin for this browser. You can click here to download the reference manual.

install.packages("covafillr")

0.4.2 by Christoffer Moesgaard Albertsen, 3 months ago


https://github.com/calbertsen/covafillr


Report a bug at https://github.com/calbertsen/covafillr/issues


Browse source code at https://github.com/cran/covafillr


Authors: Christoffer Moesgaard Albertsen [aut, cre]


Documentation:   PDF Manual  


BSD_2_clause + file LICENSE license


Imports methods, ggplot2, stats

Suggests TMB, rjags, inline

Linking to RcppEigen


See at CRAN