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'.


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")

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.

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')

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

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.

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,h = 10.0,p = 5L)

where h is the bandwith and p is the polynomial degree. Information about the class can be extracted (and changed) by the following functions:

cf$getDim()
## [1] 1
cf$getDegree()
## [1] 5
cf$getBandwith()
## [1] 10
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)

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,h = 10.0,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)

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.04614941  0.01605094

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.7055
obj$fn(c(0))
## [1] -0.04614941
obj$fn(c(-1))
## [1] -0.04697095
obj$gr()
## outer mgc:  3.593036

##           [,1]
## [1,] -3.593036
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.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.1 by Christoffer Moesgaard Albertsen, 5 months ago


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