spBayes: Interfaces for spBayes package for data science pipelines.

Description Usage Arguments Details Value Author(s) Examples

Description

Interfaces to spBayes functions that can be used in a pipeline implemented by magrittr.

Usage

1
2
3

Arguments

data

data frame, tibble, list, ...

...

Other arguments passed to the corresponding interfaced function.

Details

Interfaces call their corresponding interfaced function.

Value

Object returned by interfaced function.

Author(s)

Roberto Bertolusso

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
## Not run: 
library(intubate)
library(magrittr)
library(spBayes)

## NOTE: Just interfacing 3 functions as a proof of concept. However,
##       there are so many things declared and used, that I am not sure if
##       everything could be passed on a pipeline (perhaps with intuBags down the line).
##       I may get back to it later.

## ntbt_bayesGeostatExact: Simple Bayesian spatial linear model with fixed semivariogram parameters
data(FORMGMT.dat)

n <- nrow(FORMGMT.dat)
p <- 5 ##an intercept an four covariates

n.samples <- 50

phi <- 0.0012

coords <- cbind(FORMGMT.dat$Longi, FORMGMT.dat$Lat)
coords <- coords*(pi/180)*6378

beta.prior.mean <- rep(0, times=p)
beta.prior.precision <- matrix(0, nrow=p, ncol=p)

alpha <- 1/1.5

sigma.sq.prior.shape <- 2.0
sigma.sq.prior.rate <- 10.0


## Original function to interface
bayesGeostatExact(Y ~ X1 + X2 + X3 + X4, data = FORMGMT.dat,
                  n.samples = n.samples,
                  beta.prior.mean = beta.prior.mean,
                  beta.prior.precision = beta.prior.precision,
                  coords = coords, phi = phi, alpha = alpha,
                  sigma.sq.prior.shape = sigma.sq.prior.shape,
                  sigma.sq.prior.rate = sigma.sq.prior.rate)

## The interface puts data as first parameter
ntbt_bayesGeostatExact(FORMGMT.dat, Y ~ X1 + X2 + X3 + X4,
                       n.samples = n.samples,
                       beta.prior.mean = beta.prior.mean,
                       beta.prior.precision = beta.prior.precision,
                       coords = coords, phi = phi, alpha = alpha,
                       sigma.sq.prior.shape = sigma.sq.prior.shape,
                       sigma.sq.prior.rate = sigma.sq.prior.rate)

## so it can be used easily in a pipeline.
FORMGMT.dat %>%
  ntbt_bayesGeostatExact(Y ~ X1 + X2 + X3 + X4,
                         n.samples = n.samples,
                         beta.prior.mean = beta.prior.mean,
                         beta.prior.precision = beta.prior.precision,
                         coords = coords, phi = phi, alpha = alpha,
                         sigma.sq.prior.shape = sigma.sq.prior.shape,
                         sigma.sq.prior.rate = sigma.sq.prior.rate)




## ntbt_bayesLMConjugate: Simple Bayesian linear model via the Normal/inverse-Gamma conjugate
n <- nrow(FORMGMT.dat)
p <- 7 ##an intercept and six covariates

n.samples <- 500

## Below we demonstrate the conjugate function in the special case
## with improper priors. The results are the same as for the above,
## up to MC error. 
beta.prior.mean <- rep(0, times=p)
beta.prior.precision <- matrix(0, nrow=p, ncol=p)

prior.shape <- -p/2
prior.rate <- 0

## Original function to interface
bayesLMConjugate(Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = FORMGMT.dat,
                 n.samples, beta.prior.mean,
                 beta.prior.precision,
                 prior.shape, prior.rate)

## The interface puts data as first parameter
ntbt_bayesLMConjugate(FORMGMT.dat, Y ~ X1 + X2 + X3 + X4 + X5 + X6,
                      n.samples, beta.prior.mean,
                      beta.prior.precision,
                      prior.shape, prior.rate)

## so it can be used easily in a pipeline.
FORMGMT.dat %>%
  ntbt_bayesLMConjugate(Y ~ X1 + X2 + X3 + X4 + X5 + X6,
                        n.samples, beta.prior.mean,
                        beta.prior.precision,
                        prior.shape, prior.rate)


## ntbt_spDynLM: Function for fitting univariate Bayesian dynamic
##               space-time regression models
data("NETemp.dat")
ne.temp <- NETemp.dat

set.seed(1)

##take a chunk of New England
ne.temp <- ne.temp[ne.temp[,"UTMX"] > 5500000 & ne.temp[,"UTMY"] > 3000000,]

##subset first 2 years (Jan 2000 - Dec. 2002)
y.t <- ne.temp[,4:27]
N.t <- ncol(y.t) ##number of months
n <- nrow(y.t) ##number of observation per months

##add some missing observations to illistrate prediction
miss <- sample(1:N.t, 10)
holdout.station.id <- 5
y.t.holdout <- y.t[holdout.station.id, miss]
y.t[holdout.station.id, miss] <- NA

##scale to km
coords <- as.matrix(ne.temp[,c("UTMX", "UTMY")]/1000)
max.d <- max(iDist(coords))

##set starting and priors
p <- 2 #number of regression parameters in each month

starting <- list("beta"=rep(0,N.t*p), "phi"=rep(3/(0.5*max.d), N.t),
                 "sigma.sq"=rep(2,N.t), "tau.sq"=rep(1, N.t),
                 "sigma.eta"=diag(rep(0.01, p)))

tuning <- list("phi"=rep(5, N.t)) 

priors <- list("beta.0.Norm"=list(rep(0,p), diag(1000,p)),
               "phi.Unif"=list(rep(3/(0.9*max.d), N.t), rep(3/(0.05*max.d), N.t)),
               "sigma.sq.IG"=list(rep(2,N.t), rep(10,N.t)),
               "tau.sq.IG"=list(rep(2,N.t), rep(5,N.t)),
               "sigma.eta.IW"=list(2, diag(0.001,p)))

##make symbolic model formula statement for each month
mods <- lapply(paste(colnames(y.t),'elev',sep='~'), as.formula)

n.samples <- 10 # in original example it is 2000.

## Original function to interface
spDynLM(mods, data=cbind(y.t,ne.temp[,"elev",drop=FALSE]), coords=coords,
        starting=starting, tuning=tuning, priors=priors, get.fitted =TRUE,
        cov.model="exponential", n.samples=n.samples, n.report=25) 

## The interface puts data as first parameter
ntbt_spDynLM(cbind(y.t,ne.temp[,"elev",drop=FALSE]), mods, coords=coords,
             starting=starting, tuning=tuning, priors=priors, get.fitted =TRUE,
             cov.model="exponential", n.samples=n.samples, n.report=25) 

## so it can be used easily in a pipeline.
cbind(y.t,ne.temp[,"elev",drop=FALSE]) %>%
  ntbt_spDynLM(mods, coords=coords,
               starting=starting, tuning=tuning, priors=priors, get.fitted =TRUE,
               cov.model="exponential", n.samples=n.samples, n.report=25) 

## End(Not run)

intubate documentation built on May 2, 2019, 2:46 p.m.