binaryGP_fit {binaryGP}  R Documentation 
Binary Gaussian Process (with/without timeseries)
Description
The function fits Gaussian process models with binary response. The models can also include timeseries term if a timeseries binary response is observed. The estimation methods are based on PQL/PQPL and REML (for mean function and correlation parameter, respectively).
Usage
binaryGP_fit(X, Y, R = 0, L = 0, corr = list(type = "exponential", power =
2), nugget = 1e10, constantMean = FALSE, orthogonalGP = FALSE,
converge.tol = 1e10, maxit = 15, maxit.PQPL = 20, maxit.REML = 100,
xtol_rel = 1e10, standardize = FALSE, verbose = TRUE)
Arguments
X 
a design matrix with dimension 
Y 
a response matrix with dimension 
R 
a positive integer specifying the order of autoregression only if timeseries is included. The default is 1. 
L 
a positive integer specifying the order of interaction between 
corr 
a list of parameters for the specifing the correlation to be used. Either exponential correlation function or Matern correlation function can be used. See 
nugget 
a positive value to use for the nugget. If 
constantMean 
logical. 
orthogonalGP 
logical. 
converge.tol 
convergence tolerance. It converges when relative difference with respect to 
maxit 
a positive integer specifying the maximum number of iterations for estimation to be performed before the estimates are convergent. The default is 15. 
maxit.PQPL 
a positive integer specifying the maximum number of iterations for PQL/PQPL estimation to be performed before the estimates are convergent. The default is 20. 
maxit.REML 
a positive integer specifying the maximum number of iterations in 
xtol_rel 
a postive value specifying the convergence tolerance for 
standardize 
logical. If 
verbose 
logical. If 
Details
Consider the model
logit(p_t(x))=\eta_t(x)=\sum^R_{r=1}\varphi_ry_{tr}(\mathbf{x})+\alpha_0+\mathbf{x}'\boldsymbol{\alpha}+\sum^L_{l=1}\boldsymbol{\gamma}_l\textbf{x}y_{tl}(\mathbf{x})+Z_t(\mathbf{x}),
where p_t(x)=Pr(y_t(x)=1)
and Z_t(\cdot)
is a Gaussian process with zero mean, variance \sigma^2
, and correlation function R_{\boldsymbol{\theta}}(\cdot,\cdot)
with unknown correlation parameters \boldsymbol{\theta}
. The power exponential correlation function is used and defined by
R_{\boldsymbol{\theta}}(\mathbf{x}_i,\mathbf{x}_j)=\exp\{\sum^{d}_{l=1}\frac{(x_{il}x_{jl})^p}{\theta_l} \},
where p
is given by power
. The parameters in the mean functions including \varphi_r,\alpha_0,\boldsymbol{\alpha},\boldsymbol{\gamma}_l
are estimated by PQL/PQPL, depending on whether the mean function has the timeseries structure. The parameters in the Gaussian process including \boldsymbol{\theta}
and \sigma^2
are estimated by REML. If orthogonalGP
is TRUE
, then orthogonal covariance function (Plumlee and Joseph, 2016) is employed. More details can be seen in Sung et al. (unpublished).
Value
alpha_hat 
a vector of coefficient estimates for the linear relationship with X. 
varphi_hat 
a vector of autoregression coefficient estimates. 
gamma_hat 
a vector of the interaction effect estimates. 
theta_hat 
a vector of correlation parameters. 
sigma_hat 
a value of sigma estimate (standard deviation). 
nugget_hat 
if 
orthogonalGP 

X 
data 
Y 
data 
R 
order of autoregression. 
L 
order of interaction between X and previous Y. 
corr 
a list of parameters for the specifing the correlation to be used. 
Model.mat 
model matrix. 
eta_hat 
eta_hat. 
standardize 

mean.x 
mean of each column of 
scale.x 

Author(s)
ChihLi Sung <iamdfchile@gmail.com>
See Also
predict.binaryGP
for prediction of the binary Gaussian process, print.binaryGP
for the list of estimation results, and summary.binaryGP
for summary of significance results.
Examples
library(binaryGP)
##### Testing function: cos(x1 + x2) * exp(x1*x2) with TT sequences #####
##### Thanks to Sonja Surjanovic and Derek Bingham, Simon Fraser University #####
test_function < function(X, TT)
{
x1 < X[,1]
x2 < X[,2]
eta_1 < cos(x1 + x2) * exp(x1*x2)
p_1 < exp(eta_1)/(1+exp(eta_1))
y_1 < rep(NA, length(p_1))
for(i in 1:length(p_1)) y_1[i] < rbinom(1,1,p_1[i])
Y < y_1
P < p_1
if(TT > 1){
for(tt in 2:TT){
eta_2 < 0.3 * y_1 + eta_1
p_2 < exp(eta_2)/(1+exp(eta_2))
y_2 < rep(NA, length(p_2))
for(i in 1:length(p_2)) y_2[i] < rbinom(1,1,p_2[i])
Y < cbind(Y, y_2)
P < cbind(P, p_2)
y_1 < y_2
}
}
return(list(Y = Y, P = P))
}
set.seed(1)
n < 30
n.test < 10
d < 2
X < matrix(runif(d * n), ncol = d)
##### without timeseries #####
Y < test_function(X, 1)$Y ## Y is a vector
binaryGP.model < binaryGP_fit(X = X, Y = Y)
print(binaryGP.model) # print estimation results
summary(binaryGP.model) # significance results
##### with timeseries, lag 1 #####
Y < test_function(X, 10)$Y ## Y is a matrix with 10 columns
binaryGP.model < binaryGP_fit(X = X, Y = Y, R = 1)
print(binaryGP.model) # print estimation results
summary(binaryGP.model) # significance results