A generalized Gaussian process model for computer experiments with binary time series

TL;DR

This paper introduces a generalized Gaussian process model tailored for binary time series data in computer experiments, enabling better analysis of complex biological processes like cell adhesion.

Contribution

The paper develops a flexible GP model for binary responses with a novel mean function for time series, along with asymptotic theory and predictive methods.

Findings

Model effectively captures biological information in cell adhesion data.

Simulation studies demonstrate the model's predictive accuracy.

Application reveals insights not observable in lab experiments.

Abstract

Non-Gaussian observations such as binary responses are common in some computer experiments. Motivated by the analysis of a class of cell adhesion experiments, we introduce a generalized Gaussian process model for binary responses, which shares some common features with standard GP models. In addition, the proposed model incorporates a flexible mean function that can capture different types of time series structures. Asymptotic properties of the estimators are derived, and an optimal predictor as well as its predictive distribution are constructed. Their performance is examined via two simulation studies. The methodology is applied to study computer simulations for cell adhesion experiments. The fitted model reveals important biological information in repeated cell bindings, which is not directly observable in lab experiments.

Tables

Table 1. Table 1 : Estimation of linear coefficients. The values are the average estimates over 100 replicates, while the values in parentheses are the standard deviation of the estimates. The parameter settings are α 0 = 0.5 , α 1 = − 3 , α 2 = 2 , α 3 = − 2 , α 4 = 1 , α 5 = 0.5 formulae-sequence subscript 𝛼 0 0.5 formulae-sequence subscript 𝛼 1 3 formulae-sequence subscript 𝛼 2 2 formulae-sequence subscript 𝛼 3 2 formulae-sequence subscript 𝛼 4 1 subscript 𝛼 5 0.5 \alpha_{0}=0.5,\alpha_{1}=-3,\alpha_{2}=2,\alpha_{3}=-2,\alpha_{4}=1,\alpha_{5}=0.5 , and φ 1 = 0.8 subscript 𝜑 1 0.8 \varphi_{1}=0.8 .

$n$	$T$	Corr	${\hat{\alpha }}_0$	${\hat{\alpha }}_1$	${\hat{\alpha }}_2$	${\hat{\alpha }}_3$	${\hat{\alpha }}_4$	${\hat{\alpha }}_5$	${\hat{\phi }}_1$
200	20	PE	0.46	$-2.71$	1.82	$-1.82$	0.91	0.46	0.72
			(0.11)	(0.15)	(0.13)	(0.12)	(0.09)	(0.10)	(0.11)
		OGP	0.48	$-2.77$	1.84	$-1.85$	0.91	0.46	0.71
			(0.12)	(0.12)	(0.10)	(0.09)	(0.08)	(0.08)	(0.10)
		IMP (%)	4	2	1	1.5	0	0	$-1.25$
200	50	PE	0.45	$-2.68$	1.80	$-1.79$	0.90	0.46	0.70
			(0.07)	(0.10)	(0.09)	(0.08)	(0.07)	(0.06)	(0.07)
		OGP	0.47	$-2.75$	1.83	$-1.83$	0.92	0.46	0.71
			(0.08)	(0.09)	(0.06)	(0.06)	(0.06)	(0.05)	(0.08)
		IMP (%)	4	2.3	1.5	2	2	0	1.25
500	20	PE	0.47	$-2.76$	1.82	$-1.83$	0.91	0.48	0.73
			(0.09)	(0.14)	(0.12)	(0.10)	(0.08)	(0.07)	(0.07)
		OGP	0.49	$-2.81$	1.89	$-1.88$	0.95	0.46	0.74
			(0.08)	(0.09)	(0.07)	(0.06)	(0.06)	(0.06)	(0.07)
		IMP (%)	4	1.7	3.5	2.5	4	4	1.25
500	50	PE	0.45	$-2.75$	1.83	$-1.83$	0.92	0.45	0.74
			(0.06)	(0.07)	(0.06)	(0.06)	(0.06)	(0.05)	(0.04)
		OGP	0.47	$-2.80$	1.87	$-1.87$	0.93	0.47	0.75
			(0.06)	(0.05)	(0.04)	(0.04)	(0.03)	(0.03)	(0.04)
		IMP (%)	4	1.7	2	2	1	2	1.25

Table 2. Table 2 : Estimation of correlation parameters and variance. The values are the average estimates over 100 replicates, while the values in parentheses are the standard deviation of the estimates. The parameter settings are θ 1 = 0.5 , θ 2 = 1.0 , θ 3 = 1.5 , θ 4 = 2 , θ 5 = 2.5 formulae-sequence subscript 𝜃 1 0.5 formulae-sequence subscript 𝜃 2 1.0 formulae-sequence subscript 𝜃 3 1.5 formulae-sequence subscript 𝜃 4 2 subscript 𝜃 5 2.5 \theta_{1}=0.5,\theta_{2}=1.0,\theta_{3}=1.5,\theta_{4}=2,\theta_{5}=2.5 , and σ 2 = 1 superscript 𝜎 2 1 \sigma^{2}=1 .

$n$	$T$	${\hat{\theta }}_1$	${\hat{\theta }}_2$	${\hat{\theta }}_3$	${\hat{\theta }}_4$	${\hat{\theta }}_5$	${\hat{\sigma }}^2$
200	20	0.86	1.80	2.35	3.30	4.10	0.82
		(0.81)	(1.13)	(1.41)	(1.76)	(1.85)	(0.07)
200	50	0.65	1.55	2.38	3.01	3.80	0.79
		(0.16)	(0.63)	(1.12)	(1.25)	(1.49)	(0.05)
500	20	0.61	1.17	1.93	2.66	3.24	0.87
		(0.16)	(0.25)	(0.54)	(0.96)	(1.17)	(0.05)
500	50	0.57	1.16	1.78	2.37	3.11	0.87
		(0.08)	(0.16)	(0.35)	(0.39)	(0.68)	(0.03)

Table 3. Table 3 : The comparison of RMSPEs.

	$n=200$	$n=200$	$n=500$	$n=500$
	$T=20$	$T=50$	$T=20$	$T=50$
RMSPE	0.1188	0.1193	0.1058	0.1053
	(0.0066)	(0.0058)	(0.0060)	(0.0045)

Table 4. Table 4 : Control variables in cell adhesion frequency assay experiments.

Variable	Description	Range
$x_{K_{f,p}}$	on-rate enhancement of activated TCRs	(1,100)
$x_{K_{r,p}}$	off-rate enhancement of activated TCRs	(0.1,100)
$x_{T_{half}}$	half-life of cluster activation	(0.1,10)
$x_{Tc}$	cell-cell contact time	(0.1,10)
$x_{Tw}$	waiting time in between contacts	(0.1,10)
$x_{Kc}$	kinetic proofreading modification rate for activation of cluster	(0.1,10)

Table 5. Table 5 : Estimation results.

		Standard
	Value	deviation	$Z$ score	p value
${\hat{\phi }}_1$	0.14	0.02	7.96	0.0000
${\hat{\alpha }}_1$	0.13	0.02	5.96	0.0000
${\hat{\alpha }}_2$	0.37	0.02	17.82	0.0000
${\hat{\alpha }}_3$	0.00	0.02	0.08	0.9331
${\hat{\alpha }}_4$	0.47	0.02	22.54	0.0000
${\hat{\alpha }}_5$	-0.08	0.02	-3.79	0.0001
${\hat{\alpha }}_6$	0.15	0.02	6.86	0.0000
${\hat{\gamma }}_1$	0.16	0.03	5.2	0.0000
${\hat{\gamma }}_2$	0.23	0.03	7.68	0.0000
${\hat{\gamma }}_3$	-0.09	0.03	-2.92	0.0035
${\hat{\gamma }}_4$	0.23	0.03	7.28	0.0000
${\hat{\gamma }}_5$	-0.17	0.03	-5.47	0.0000
${\hat{\gamma }}_6$	0.36	0.03	11.83	0.0000

Table 6. Table 6 : Comparison of prediction performance using proper scoring rules. The values are the medians over the 10 folds. Larger value indicates better prediction. binaryGP: proposed method, glm: logistic regression, glm_ts: logistic regression with time-series mean function, kernlab: Bayesian generalized GP, and GPFDA: functional Gaussian process model.

Scoring rule	binaryGP	glm	glm_ts	kernlab	GPFDA
Brier	-0.2092	-0.2193	-0.2199	-0.2197	-0.2145
Spherical	0.7632	0.7499	0.7490	0.7491	0.7545
Logarithmic	-0.6099	-0.6310	-0.6323	-0.6323	-0.6166
Zero-one ($c=1/2$)	0.3688	0.3517	0.3542	0.3246	0.3650

Equations

$$\text{logit}(p(\mathbf{x}))={\alpha }_0+{\mathbf{x}}^{\prime }\mathbf{\alpha }+Z(\mathbf{x}),$$

$$R_{\mathbf{\theta }}({\mathbf{x}}_i,{\mathbf{x}}_j)=\exp \left\{-\sum _{l=1}^d\frac{(x_{il}-x_{jl}{)}^p}{{\theta }_l}\right\},$$

$$\text{logit}(p_t(\mathbf{x}))={\eta }_t(\mathbf{x})=\sum _{r=1}^R{\phi }_r{y}_{t-r}(\mathbf{x})+{\alpha }_0+{\mathbf{x}}^{\prime }\mathbf{\alpha }+\sum _{l=1}^L{\mathbf{\gamma }}_l\text{x}{y}_{t-l}(\mathbf{x})+Z_t(\mathbf{x}),$$

$$\text{logit}(\mathbf{p})=\mathbf{X}\mathbf{\beta }+\mathbf{Z},\mathbf{Z}\sim \mathcal{N}(0_N,\Sigma (\mathbf{\omega })),$$

$$\Sigma (\mathbf{\omega })={\sigma }^2{\mathbf{R}}_{\mathbf{\theta }}\otimes I_T$$

$$PL(\mathbf{\beta }|\mathbf{Z})=\prod _{i=1}^n\prod _{t=1}^T(p_{it}(\mathbf{\beta }|\mathbf{Z}){)}^{y_{it}}(1-p_{it}(\mathbf{\beta }|\mathbf{Z}){)}^{1-y_{it}},$$

$$|\Sigma (\mathbf{\omega }){|}^{-1/2}\int \exp \{\log PL(\mathbf{\beta }|\mathbf{Z})-\frac{1}{2}{\mathbf{Z}}^{\prime }\Sigma (\mathbf{\omega }{)}^{-1}\mathbf{Z}\}d\mathbf{Z}.$$

$$L(\mathbf{\beta }|\mathbf{Z})=\prod _{i=1}^n(p_{i1}(\mathbf{\beta }|\mathbf{Z}){)}^{y_{i1}}(1-p_{i1}(\mathbf{\beta }|\mathbf{Z}){)}^{1-y_{i1}}$$

$$|\Sigma (\mathbf{\omega }){|}^{-1/2}\int \exp \{\log L(\mathbf{\beta }|\mathbf{Z})-\frac{1}{2}{\mathbf{Z}}^{\prime }\Sigma (\mathbf{\omega }{)}^{-1}\mathbf{Z}\}d\mathbf{Z},$$

$$-\frac{1}{2}\log |I_n$$

$$\sum _{i=1}^n\sum _{t=1}^T\left(y_{it}\log \frac{p_{it}(\mathbf{\beta }|\tilde{\mathbf{Z}})}{1-p_{it}(\mathbf{\beta }|\tilde{\mathbf{Z}})}+\log (1-p_{it}(\mathbf{\beta }|\tilde{\mathbf{Z}}))\right)-\frac{1}{2{\sigma }^2}{\tilde{\mathbf{Z}}}^{\prime }({\mathbf{R}}_{\mathbf{\theta }}\otimes I_T{)}^{-1}\tilde{\mathbf{Z}},$$

$$\sum _{i=1}^n\sum _{t=1}^T{X}_{it}(y_{it}-p_{it}(\mathbf{\beta },\mathbf{\omega }))=0$$

$$\sum _{i=1}^n\sum _{t=1}^T{\mathbf{e}}_{it}(y_{it}-p_{it}(\mathbf{\beta },\mathbf{\omega }))=({\mathbf{R}}_{\mathbf{\theta }}\otimes I_T{)}^{-1}\mathbf{Z}/{\sigma }^2,$$

$$({\mathbf{X}}^{\prime }\mathbf{V}(\mathbf{\omega }{)}^{-1}\mathbf{X})\mathbf{\beta }={\mathbf{X}}^{\prime }\mathbf{V}(\mathbf{\omega }{)}^{-1}\tilde{\mathbf{\eta }},$$

$$\hat{\mathbf{Z}}={\sigma }^2({\mathbf{R}}_{\mathbf{\theta }}\otimes I_T)\mathbf{V}(\mathbf{\omega }{)}^{-1}(\tilde{\mathbf{\eta }}-{\mathbf{X}}^{\prime }\hat{\mathbf{\beta }})$$

$$L(\mathbf{\omega })=\frac{N-m}{2}\log (2\pi )-\frac{1}{2}\log (|{\mathbf{X}}^{\prime }\mathbf{X}|)+\frac{1}{2}\log (|\mathbf{V}(\mathbf{\omega })|)+\frac{1}{2}\log (|{\mathbf{X}}^{\prime }\mathbf{V}(\mathbf{\omega }{)}^{-1}\mathbf{X}|)+\frac{1}{2}{\tilde{\mathbf{\eta }}}^{\prime }\Pi (\mathbf{\omega })\tilde{\mathbf{\eta }},$$

$$\sqrt{N}(\hat{\mathbf{\beta }}-\mathbf{\beta })={\Lambda }_N^{-1}\frac{1}{\sqrt{N}}{S}_N(\mathbf{\beta },\mathbf{\omega })+o_p(1)$$

$$\sqrt{N}{\Lambda }_N^{1/2}(\hat{\mathbf{\beta }}-\mathbf{\beta })\stackrel{d}{⟶}\mathcal{N}(0,I_m),$$

$${\Lambda }_N=\frac{1}{N}\sum _{i=1}^n\sum _{t=1}^T{X}_{it}{X}_{it}^{\prime }{p}_{it}(\mathbf{\beta },\mathbf{\omega })(1-p_{it}(\mathbf{\beta },\mathbf{\omega })),$$

$$J_N(\hat{\mathbf{\omega }})(\hat{\mathbf{\omega }}-\mathbf{\omega })\stackrel{d}{⟶}\mathcal{N}(0,I_{d+1}).$$

$$p_s({\mathbf{x}}_{n+1})|D_{n+1,s}\sim Logitnormal(m(D_{n+1,s}),v(D_{n+1,s})),$$

$$m(D_{n+1,s})=\sum _{r=1}^R{\phi }_r{y}_{n+1,s-r}+{\alpha }_0+{\mathbf{x}}_{n+1}^{\prime }\mathbf{\alpha }+\sum _{l=1}^L{\mathbf{\gamma }}_l{\mathbf{x}}_{n+1}{y}_{n+1,s-l}+{\mathbf{r}}_{\mathbf{\theta }}^{\prime }{\mathbf{R}}_{\mathbf{\theta }}^{-1}\left(\log \frac{{\mathbf{p}}_s}{1_n-{\mathbf{p}}_s}-{\mathbf{\mu }}_s\right),$$

$$p({\mathbf{x}}_{n+1})|D_{n+1}\sim Logitnormal(m(D_{n+1}),v(D_{n+1})),$$

$$\hat{p}=\hat{p}(D)=\arg \underset{\eta }{\min }{\mathbb{E}}_F[(p-\eta {)}^2],$$

$$\mathbb{E}\left[p_s({\mathbf{x}}_{n+1})|D_{n+1,s}\right]$$

$$\text{with variance}\mathbb{V}\left[p_s({\mathbf{x}}_{n+1})|D_{n+1,s}\right]$$

$$\frac{\exp \{m(D_{n+1,s})+z_q\sqrt{v(D_{n+1,s})}\}}{1+\exp \{m(D_{n+1,s})+z_q\sqrt{v(D_{n+1,s})}\}},$$

$${\hat{p}}_s({\mathbf{x}}_{n+1})=\mathbb{E}\left[p_s({\mathbf{x}}_{n+1})|\mathbf{Y}\right]={\mathbb{E}}_{\mathbf{p}|\mathbf{Y}}\left[\kappa (m(D_{n+1,s}),v(D_{n+1,s}))|\mathbf{Y}\right]$$

$$\mathbb{V}\left[p_s({\mathbf{x}}_{n+1})|\mathbf{Y}\right]={\mathbb{E}}_{\mathbf{p}|\mathbf{Y}}\left[\tau (m(D_{n+1,s}),v(D_{n+1,s}))|\mathbf{Y}\right]+{\mathbb{V}}_{\mathbf{p}|\mathbf{Y}}\left[\kappa (m(D_{n+1,s}),v(D_{n+1,s}))|\mathbf{Y}\right],$$

$${\mathbb{E}}_{\mathbf{p}|\mathbf{Y}}\left[\kappa (m(D_{n+1,s}),v(D_{n+1,s}))|\mathbf{Y}\right]\approx \frac{1}{J}\sum _{j=1}^J\kappa (m(D_{n+1,s}^{(j)}),v(D_{n+1,s}^{(j)})),$$

$$\text{RMSPE}={\left(\frac{1}{n_{\text{test}}T}\sum _{i=1}^{n_{\text{test}}}\sum _{t=1}^T(p_t({\mathbf{x}}_i^{\ast })-{\hat{p}}_t({\mathbf{x}}_i^{\ast }){)}^2\right)}^{1/2},$$

Full text

Keywords: Computer experiment, Gaussian process model, Single molecule experiment, Uncertainty quantification

1 Introduction

Cell adhesion plays an important role in many physiological and pathological processes. This research is motivated by the analysis of a class of cell adhesion experiments called micropipette adhesion frequency assays, which is a method for measuring the kinetic rates between molecules in their native membrane environment. In a micropipette adhesion frequency assay, a red blood coated in a specific ligand is brought into contact with cell containing the native receptor for a predetermined duration, then retracted. The output of interest is binary, indicating whether a controlled contact results in adhesion. If there is an adhesion between molecules at the end of contact, retraction will stretch the red cell. If no adhesion resulted, the red cell will not be stretched. The kinetics of the molecular interaction can be derived through many repeated trials. In theory, these contacts should be independent Bernoulli trials. However, there is a memory effect in the repeated tests and the quantification of such a memory effect is scientifically important (Zarnitsyna et al.,, 2007; Hung et al.,, 2008).

A cost-effective way to study the repeated adhesion frequency assays is through computer experiments, which study real systems using complex mathematical models and numerical tools such as finite element analysis (Santner et al.,, 2003). They have been widely used as alternatives to physical experiments or observations, especially for the study of complex systems. For cell adhesion, performing physical experiments (i.e., lab work) is time-consuming and often involves complicated experimental manipulation. Therefore, instead of performing the experiments only based on the actual lab work, computer simulations based on mathematical models are conducted to provide an efficient way to examine the complex mechanisms behind the adhesion.

The analysis of computer experiments has three objectives: (i) to build a model that captures the nonlinear relationship between inputs and outputs; (ii) to estimate the unknown parameters in the model and deduce properties of the estimators; (iii) to provide an optimal predictor for untried input settings, also called “emulator” or “surrogate model”, and quantify its predictive uncertainty (Sacks et al.,, 1989; Santner et al.,, 2003). This objective (iii) is crucial because computer simulations are generally expensive or time-consuming to perform and therefore the emulators based on computer simulations are used as surrogates to perform sensitivity analysis, process optimization, calibration, etc. In particular, it is critical for calibration problems in which the emulators and physical experiments are integrated so that some unknown calibration parameters can be estimated. In the literature, Gaussian process (GP) model, use of which achieves the three objectives, is widely used for the analysis of computer experiments. A GP model accommodates nonlinearity using GP and provides an optimal predictor with an interpolation property. The applications of GP can be found in many fields in science and engineering.

The conventional GP models are developed for continuous outputs with a Gaussian assumption, which does not hold in some scientific studies. For example, the focus of the cell adhesion frequency assays is to elicit the relationship between the setting of kinetic parameters/covariates and the adhesion score, which is binary. For binary outputs, the Gaussian assumption is not valid and GP models cannot be directly applied. Binary outputs are common in computer experiments, but the extensions of GP models to non-Gaussian cases have received scant attention in computer experiment literature. Although there are intensive studies of generalized GP models for non-Gaussian data in machine learning and spatial statistics literature, such as Williams and Barber, (1998), Zhang, (2002), Rasmussen and Williams, (2006), Nickisch and Rasmussen, (2008) and Wang and Shi, (2014), the asymptotic properties of estimators have not been systematically studied. Moreover, an analogy to the GP predictive distribution for binary data is important for uncertainty quantification in computer experiments, which has not yet been developed to the best of our knowledge.

Apart from the non-Gaussian responses, analysis of the repeated cell adhesion frequency assays poses another challenge, namely, how to incorporate a time series structure with complex interaction effects. It was discovered that cells appear to have the ability to remember the previous adhesion events and such a memory has an impact on the future adhesion behaviors (Zarnitsyna et al.,, 2007; Hung et al.,, 2008). The quantification of the memory effect and how it interacts with the settings of the kinetic parameters in the binary time series are important but cannot be obtained by direct application of the conventional GP models. To consider the time series structure, a common practice is to construct a spatial-temporal model. However, a separable correlation function (e.g., Gelfand et al., (2004); Conti and O’Hagan, (2010)) in which space and time are assumed to be independent is often implemented as a convenient way to address the computational issue. As a result, the estimation of interaction between space and time, which is of major interest here, is not allowed for. Even in the cases where nonseparable correlation functions (e.g., Gelfand et al., (2004); Fricker et al., (2013)) are implemented, the interaction effect is still not easily interpretable. Therefore, a new model that can model binary time series and capture interaction effects is called for.

To achieve the objectives in the analysis of computer experiments and overcome the aforementioned limitations with binary time series outputs, we introduce a new class of models in this article. The idea is to generalize GP models to non-Gaussian responses and incorporate a flexible mean function that can estimate the time series structure and its interaction with the input variables. In particular, we focus on binary responses and introduce a new model which is analogous to the GP model with an optimal interpolating predictor. Rigorous studies of estimation, prediction, and inference are required for the proposed model and the derivations are complicated by the nature of binary responses and the dependency of time series. Since binary responses with serial correlations can be observed in computer experiments, the proposed method can be readily applicable to other fields beyond cell biology. For example, in manufacturing industry computer simulations are often conducted for the failure analysis where the outputs of interest are binary, i.e., failure or success (Yan et al.,, 2009). Examples can also be found in other biological problems where binary outputs are observed and evolve in time, such as neuron firing simulations, cell signaling pathways, gene transcription, and recurring diseases (Gerstner et al.,, 1997; Mayrhofer et al.,, 2002). The proposed method can also be broadly applied beyond computer experiments. In many scientific experiments, such as medical research and social studies, binary repeated measurements are commonly observed with serial correlations. In these situations, the proposed method can be implemented to provide a flexible nonlinear model that quantifies the correlation structure and explains the complex relationship between inputs and binary outputs. More examples can be found in functional data analysis, longitudinal data analysis, and machine learning.

The remainder of this article is organized as follows. The new class of models is discussed in Section 2. In Section 3 and 4, asymptotic properties of the estimators are derived and the predictive distributions are constructed. Finite sample performance is demonstrated by simulations in Section 5. In Section 6, the proposed method is illustrated with the analysis of computer experiments for cell adhesion frequency assays. Concluding remarks are given in Section 7. Mathematical proofs and algorithms are provided in the online supplement. An implementation for our method can be found in binaryGP (Sung,, 2017) in R (R Core Team,, 2015).

2 Model

2.1 Generalized Gaussian process models for binary response

We first introduce a model for binary responses in computer experiments which is analogous to the conventional GP models for continuous outputs. Suppose a computer experiment has a $d$-dimensional input setting $\mathbf{x}=(x_{1},\ldots,x_{d})^{\prime}$ and for each setting the binary output is denoted by $y(\mathbf{x})$ and randomly generated from a Bernoulli distribution with probability $p(\mathbf{x})$. Using a logistic link function, the Gaussian process model for binary data can be written as

[TABLE]

where $p(\mathbf{x})=\mathbb{E}[y(\mathbf{x})]$, $\alpha_{0}$ and $\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{d})^{\prime}$ are the intercept and linear effects of the mean function of $p(\mathbf{x})$, and $Z(\cdot)$ is a zero mean Gaussian process with variance $\sigma^{2}$, correlation function $R_{\boldsymbol{\theta}}(\cdot,\cdot)$, and unknown correlation parameters $\boldsymbol{\theta}$.

Various choices of correlation functions have been discussed in the literature. For example, the power exponential correlation function is commonly used in the analysis of computer experiments (Santner et al.,, 2003):

[TABLE]

where $\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{d})$, the power $p$ controls the smoothness of the output surface, and the parameter $\theta_{l}$ controls the decay of correlation with respect to the distance between $x_{il}$ and $x_{jl}$. Recent studies have shown that a careful selection of the correlation function, such as orthogonal Gaussian processes proposed by Plumlee and Joseph, (2018), can resolve the identifiability issue in the estimation of Gaussian process models (Hodges and Reich,, 2010; Paciorek,, 2010; Tuo and Wu,, 2015). This is particularly important in the application of calibration problems where the parameter estimation plays a significant role. Depending on the objectives of the studies, different correlation functions can be incorporated into the proposed model and the theoretical results developed herein remain valid.

Similar extensions of GP models to binary outputs have been applied in many different fields. For example, when $\mathbf{x}$ represents a two-dimensional spatial domain, (1) becomes the spatial generalized linear mixed model proposed by Zhang, (2002). In a Bayesian framework, Gaussian process priors are implemented for classification problems, such as in Williams and Barber, (1998) and Gramacy and Polson, (2011). Despite successful applications of these models, theoretical studies on the estimation and prediction properties are not available. Therefore, one focus of this paper is to provide theoretical supports for the estimation and prediction in (1).

2.2 Generalized Gaussian process models for binary time series

In this section, we introduce a new model for the analysis of computer experiments with binary time series, which is an extension of (1) that takes serial correlations between binary observations into account. Suppose for each setting of a computer experiment, a sequence of binary time series outputs $\{y_{t}(\mathbf{x})\}^{T}_{t=1}$ is randomly generated from Bernoulli distributions with probabilities $\{p_{t}(\mathbf{x})\}^{T}_{t=1}$. A generalized Gaussian process model for binary time series can be written as:

[TABLE]

where $p_{t}(\mathbf{x})=\mathbb{E}[y_{t}(\mathbf{x})|H_{t}]$ is the conditional mean given the previous information $H_{t}=\{y_{t-1}(\mathbf{x}),y_{t-2}(\mathbf{x}),\ldots\}$. In model (3), $\{\varphi_{r}\}^{R}_{r=1}$ represents an autoregressive (AR) process with order $R$ and $\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{d})^{\prime}$ represents the effects of $\mathbf{x}$. The $d$-dimensional vector $\boldsymbol{\gamma}_{l}$ represents the interaction between the input and the past outputs and provides the flexibility of modeling different time series structures with different inputs. Given that the interactions between $\mathbf{x}$ and time are captured by $\textbf{x}y_{t-l}$, $Z_{t}$ is assumed to vary independently over time to reduce modeling and computational complexity. Further extensions can be made by replacing $Z_{t}(\mathbf{x})$ with a spatio-temporal Gaussian process $Z(t,\mathbf{x})$, but the computational cost will be higher. Without the Gaussian process assumption in (3), the mean function is closely related to the Zeger-Qaqish (1988) model and its extensions in Hung et al., (2008) and Benjamin et al., (2003), all of which take into account the autoregressive predictors in logistic regression.

Model (3) extends the applications of conventional GP to binary time series generated from computer experiments. The model is intuitively appealing; however, the issues of estimation, prediction, and inference are not straightforward due to the nature of binary response and the dependency structure.

3 Inference

Since model (1) can be written as a special case of model (3) when $R=0,L=0$ and $T=1$, derivations herein are mainly based on model (3) with additional discussions given for (1) when necessary.

3.1 Estimation

Given $n$ input settings $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ in a computer experiment, denote $y_{it}\equiv y_{t}(\mathbf{x}_{i})$ as the binary output generated from input $\mathbf{x}_{i}$ at time $t$, where $\mathbf{x}_{i}\in\mathbb{R}^{d}$, $i=1,...,n$, and $t=1,...,T$. Let $N$ be the total number of the outputs, i.e., $N=nT$. In addition, at each time $t$, denote $\mathbf{y}_{t}$ as an $n$-dimensional vector $\mathbf{y}_{t}=(y_{1t},...,y_{nt})^{\prime}$ with conditional mean $\boldsymbol{p}_{t}=(p_{1t},\ldots,p_{nt})^{\prime}$, where $p_{it}=\mathbb{E}(y_{it}|H_{it})$ and $H_{it}=\{y_{i,t-1},y_{i,t-2},\ldots\}$. Based on the data, model (3) can be rewritten into matrix form as follows:

[TABLE]

where $\boldsymbol{p}=(\boldsymbol{p}^{\prime}_{1},\ldots,\boldsymbol{p}^{\prime}_{T})^{\prime},\boldsymbol{\beta}=(\varphi_{1},\ldots,\varphi_{R},\alpha_{0},\boldsymbol{\alpha}^{\prime},(\boldsymbol{\gamma}_{1}^{\prime},\ldots,\boldsymbol{\gamma}_{L}^{\prime})^{\prime})^{\prime},\boldsymbol{\omega}=(\sigma^{2},\boldsymbol{\theta})^{\prime},\\ \boldsymbol{Z}=(Z_{1}(\mathbf{x}_{1}),\ldots,Z_{1}(\mathbf{x}_{n}),\ldots,Z_{T}(\mathbf{x}_{1}),\ldots,Z_{T}(\mathbf{x}_{n}))^{\prime}$, $\boldsymbol{X}$ is the model matrix $(X^{\prime}_{1},\ldots,X^{\prime}_{T})^{\prime}$, $X_{t}$ is an $n\times(1+R+d+dL)$ matrix with $i$-th row defined by $(X_{t})_{i}=(1,y_{i,t-1},\ldots,y_{i,t-R},\mathbf{x}_{i}^{\prime},\\ \mathbf{x}^{\prime}_{i}y_{i,t-1},\ldots,\mathbf{x}^{\prime}_{i}y_{i,t-L})$, and $\Sigma(\boldsymbol{\omega})$ is an $N\times N$ covariance matrix defined by

[TABLE]

with $(\boldsymbol{R}_{\boldsymbol{\theta}})_{ij}=R_{\boldsymbol{\theta}}(\mathbf{x}_{i},\mathbf{x}_{j})$. Model (1) can also be rewritten in the same way by setting $R=0,L=0$ and $T=1$.

With the presence of time series and their interaction with the input settings in model (3), we can write down the partial likelihood (PL) function (Cox,, 1972, 1975) according to the formulation of Slud and Kedem, (1994). Given the previous information $\{H_{it}\}_{i=1,\ldots,n;t=1,\ldots,N}$, the PL for $\boldsymbol{\beta}$ can be written as

[TABLE]

where $p_{it}(\boldsymbol{\beta}|\boldsymbol{Z})=\mathbb{E}_{\boldsymbol{\beta}|\boldsymbol{Z}}[y_{it}|H_{it}]$. Then, the integrated quasi-PL function for the estimation of $(\boldsymbol{\beta},\boldsymbol{\omega})$ is given by

[TABLE]

Note that, for model (1) where no time series effect is considered, (6) and (7) should be replaced by the likelihood function

[TABLE]

and the integrated quasi-likelihood function

[TABLE]

respectively. Hereafter, we provide the framework for the integrated quasi-PL function (7), but the result can be applied to the integrated quasi-likelihood function (8) by assuming $R=0,L=0$ and $T=1$.

Because of the difficulty in computing the integrated quasi-PL function, a penalized quasi-PL (PQPL) function is used as an approximation. Similar to the procedure in Breslow and Clayton, (1993), the integrated quasi-partial log-likelihood can be approximated by Laplace’s method (Barndorff-Nielsen and Cox,, 1997). Ignoring the multiplicative constant and plugging (5) in $\Sigma(\boldsymbol{\omega})$, the approximation yields

[TABLE]

where $\boldsymbol{W}$ is an $N\times N$ diagonal matrix with diagonal elements $W_{it}=p_{it}(\boldsymbol{\beta}|\tilde{\boldsymbol{Z}})(1-p_{it}(\boldsymbol{\beta}|\tilde{\boldsymbol{Z}}))$, $p_{it}(\boldsymbol{\beta}|\tilde{\boldsymbol{Z}})=\mathbb{E}_{\boldsymbol{\beta}|\tilde{\boldsymbol{Z}}}[y_{it}|H_{it}]$, and $\tilde{\boldsymbol{Z}}=\tilde{\boldsymbol{Z}}(\boldsymbol{\beta},\boldsymbol{\omega})$ is the solution of $\sum^{n}_{i=1}\sum^{T}_{t=1}\mathbf{e}_{it}(y_{it}-p_{it}(\boldsymbol{\beta}|\boldsymbol{Z}))=(\boldsymbol{R}_{\boldsymbol{\theta}}\otimes I_{T})^{-1}\boldsymbol{Z}/\sigma^{2}$, where $\mathbf{e}_{it}$ is a unit-vector where $((t-1)n+i)$-th element is one. The estimator $\hat{\boldsymbol{\beta}}$ which maximizes the PQPL function (3.1) is called maximum quasi-PL estimator. Thus, similar to the derivations in Breslow and Clayton, (1993) for score equations of a penalized quasi-likelihood function, the score equations of the PQPL function for $\boldsymbol{\beta}$ and $\boldsymbol{\omega}$ are

[TABLE]

and

[TABLE]

where $p_{it}(\boldsymbol{\beta},\boldsymbol{\omega})=\mathbb{E}_{\boldsymbol{\beta},\boldsymbol{\omega}}[y_{it}|H_{it}]$. The solution to the score equations can be efficiently obtained by an iterated weighted least squares (IWLS) approach as follows. In each step, one first solves for $\boldsymbol{\beta}$ in

[TABLE]

where $\boldsymbol{V}(\boldsymbol{\omega})=\boldsymbol{W}^{-1}+\sigma^{2}(\boldsymbol{R}_{\boldsymbol{\theta}}\otimes I_{T})$, $\boldsymbol{W}$ is an $N\times N$ diagonal matrix with diagonal elements $W_{it}=p_{it}(\boldsymbol{\beta},\boldsymbol{\omega})(1-p_{it}(\boldsymbol{\beta},\boldsymbol{\omega}))$, and $\tilde{\eta}_{it}=\log\frac{p_{it}(\boldsymbol{\beta},\boldsymbol{\omega})}{1-p_{it}(\boldsymbol{\beta},\boldsymbol{\omega})}+\frac{y_{it}-p_{it}(\boldsymbol{\beta},\boldsymbol{\omega})}{p_{it}(\boldsymbol{\beta},\boldsymbol{\omega})(1-p_{it}(\boldsymbol{\beta},\boldsymbol{\omega}))}$, and then sets

[TABLE]

and replaces $p_{it}(\boldsymbol{\beta},\boldsymbol{\omega})$ with $p_{it}(\hat{\boldsymbol{\beta}},\boldsymbol{\omega})=\left(\frac{\exp\{\boldsymbol{X}^{\prime}\hat{\boldsymbol{\beta}}+\hat{\boldsymbol{Z}}\}}{\mathbf{1}_{N}+\exp\{\boldsymbol{X}^{\prime}\hat{\boldsymbol{\beta}}+\hat{\boldsymbol{Z}}\}}\right)_{it}$.

Estimation of the correlation parameters $\boldsymbol{\theta}$ and variance $\sigma^{2}$ is obtained by the restricted maximum likelihood (REML) approach (Patterson and Thompson,, 1971) because it is known to have smaller bias comparing with the maximum likelihood approach (Patterson and Thompson,, 1974). See also Harville, (1977) and Searle et al., (2009) for details. According to Harville, (1974, 1977), the REML estimators of $\sigma^{2}$ and $\boldsymbol{\theta}$ can be solved by minimizing the following negative log-likelihood function with respect to $\boldsymbol{\omega}$,

[TABLE]

where $m=1+R+d+dL$ and $\varPi(\boldsymbol{\omega})=\boldsymbol{V}(\boldsymbol{\omega})^{-1}-\boldsymbol{V}(\boldsymbol{\omega})^{-1}\boldsymbol{X}(\boldsymbol{X}^{\prime}\boldsymbol{V}(\boldsymbol{\omega})^{-1}\boldsymbol{X})^{-1}\boldsymbol{X}^{\prime}\boldsymbol{V}(\boldsymbol{\omega})^{-1}$.

Therefore, the estimators $\hat{\boldsymbol{\beta}}$ and $\hat{\boldsymbol{\omega}}$ $(\equiv(\hat{\sigma}^{2},\hat{\boldsymbol{\theta}})^{\prime})$ can be obtained by iteratively solving (10), (11) and minimizing (12). The explicit algorithm is given in the supplementary material S1. Note that $\boldsymbol{V}(\boldsymbol{\omega})$ is a block diagonal matrix, i.e., a square matrix having main diagonal blocks square matrices such that the off-diagonal blocks are zero matrices. Therefore the computational burden for the matrix inversion of $\boldsymbol{V}(\boldsymbol{\omega})$ can be alleviated by the fact that the inverse of a block diagonal matrix is a block diagonal matrix, composed of the inversion of each block.

3.2 Asymptotic Properties

Asymptotic results are presented here to show that the estimators $\hat{\boldsymbol{\beta}},\hat{\sigma}^{2}$ and $\hat{\boldsymbol{\theta}}$ obtained in Section 3.1 are asymptotically normally distributed when $N(=nT)$ becomes sufficiently large. In the present context both $n$ and $T$ are sufficiently large. The assumptions are given in the supplementary material S2, and the proofs are stated in the supplementary material S3 and S4. These results are developed along the lines described in Hung et al., (2008) and Cressie and Lahiri, (1993, 1996).

Theorem 3.1**.**

Under assumptions S2.1 and S2.2, the maximum quasi-PL estimator for the fixed effects $\boldsymbol{\beta}$ are consistent and asymptotically normal as $N\rightarrow\infty$,

[TABLE]

and

[TABLE]

where $m$ is the size of the vector $\boldsymbol{\beta}$ (i.e., $m=1+R+d+dL$), the sample information matrix

[TABLE]

and $S_{N}(\boldsymbol{\beta},\boldsymbol{\omega})=\sum^{n}_{i=1}\sum^{T}_{t=1}X_{it}(y_{it}-p_{it}(\boldsymbol{\beta},\boldsymbol{\omega}))$.

Remark 3.2**.**

For model (1), the estimator $\hat{\boldsymbol{\beta}}$ can be obtained by minimizing the penalized quasi-likelihood (PQL) function, which can be written as (3.1) with $T=1$. Under assumption S2.1 and the application of central limit theorem, such estimator has the same asymptotic properties as in Theorem 3.1 with $N=n$.

For models (1) and (3), we have the following asymptotic properties for $\hat{\boldsymbol{\omega}}$.

Theorem 3.3**.**

Denote $[\Gamma_{N}(\boldsymbol{\omega})]_{i,j}=\partial^{2}L(\boldsymbol{\omega})/\partial\omega_{i}\partial\omega_{j}$ and $J_{N}(\boldsymbol{\omega})=[\mathbb{E}_{\boldsymbol{\omega}}\Gamma_{N}(\boldsymbol{\omega})]^{1/2}$. Then, under assumptions S2.3 and S2.4, as $N\rightarrow\infty$,

[TABLE]

Note that the asymptotic results here focus on the conditional inference of $\hat{\boldsymbol{\beta}}|\boldsymbol{\omega}$ and $\hat{\boldsymbol{\omega}}|\boldsymbol{\beta}$. Therefore, these results still hold in the presence of the orthogonal Gaussian process approach proposed by Plumlee and Joseph, (2018). Theoretically speaking, the orthogonal Gaussian process approach is expected to reduce the covariance between $\hat{\boldsymbol{\beta}}$ and $\hat{\boldsymbol{\omega}}$. However, the derivation of the joint distribution would be nontrivial. We leave this result to the future work.

4 Construction of Predictive Distribution

For computer experiments, the construction of an optimal predictor and its corresponding predictive distribution is important for uncertainty quantification, sensitivity analysis, process optimization, and calibration (Santner et al.,, 2003).

First, some notation is introduced. For some untried setting $\mathbf{x}_{n+1}$, denote the predictive probability at time $s$ by $p_{s}(\mathbf{x}_{n+1})=\mathbb{E}[y_{s}(\mathbf{x}_{n+1})|H_{s}]$, where $H_{s}=\{y_{n+1,s-1},y_{n+1,s-2},\ldots\}$. Assume that $D_{n+1,s}$ represents the “previous information” including $\{y_{n+1,s-1},y_{n+1,s-2},\ldots,\\ p_{n+1,s-1},p_{n+1,s-2},\ldots\}$ at $\mathbf{x}_{n+1}$ and $\{y_{it},p_{it}\}$, where $i=1,\ldots,n$ and $t=1,\ldots,T$. Also, let $Logitnormal(\mu,\sigma^{2})$ represent a logit-normal distribution $P$, where $P=\exp\{X\}/(1+\exp\{X\})$ and $X$ has a univariate normal distribution with $\mu$ and variance $\sigma^{2}$. Denote the first two moments of the distribution by $\mathbb{E}[P]=\kappa(\mu,\sigma^{2})$ and $\mathbb{V}[P]=\tau(\mu,\sigma^{2})$. In general, there is no closed form expression for $\kappa(\mu,\sigma^{2})$ and $\tau(\mu,\sigma^{2})$, but it can be easily computed by numerical integration such as in the package logitnorm (Wutzler,, 2012) in R (R Core Team,, 2015). More discussions on logit-normal distribution can be found in Mead, (1965); Atchison and Shen, (1980); Frederic and Lad, (2008).

We first present a lemma which shows that, given $D_{n+1,s}$, the conditional distribution of $p_{s}(\mathbf{x}_{n+1})$ in model (3) is logit-normal. This result lays the foundation for the construction of predictive distribution. The proof is given in the supplementary material S5.

Lemma 4.1**.**

For model (3), the conditional distribution of $p_{s}(\mathbf{x}_{n+1})$ can be written as

[TABLE]

where

[TABLE]

$v(D_{n+1,s})=\sigma^{2}\left(1-\boldsymbol{r}_{\boldsymbol{\theta}}^{\prime}\boldsymbol{R}_{\boldsymbol{\theta}}^{-1}\boldsymbol{r}_{\boldsymbol{\theta}}\right),$* $\boldsymbol{r}_{\boldsymbol{\theta}}=(R_{\boldsymbol{\theta}}(\mathbf{x}_{n+1},\mathbf{x}_{1}),\ldots,R_{\boldsymbol{\theta}}(\mathbf{x}_{n+1},\mathbf{x}_{n}))^{\prime},\boldsymbol{R}_{\boldsymbol{\theta}}=\{R_{\boldsymbol{\theta}}(\mathbf{x}_{i},\mathbf{x}_{j})\}$, $\boldsymbol{p}_{s}=(p_{s}(\mathbf{x}_{1}),\ldots,p_{s}(\mathbf{x}_{n}))^{\prime}$, and $(\boldsymbol{\mu}_{s})_{i}=\sum^{R}_{r=1}\varphi_{r}y_{i,s-r}+\alpha_{0}+\mathbf{x}_{i}^{\prime}\boldsymbol{\alpha}+\sum^{L}_{l=1}\boldsymbol{\gamma}_{l}\mathbf{x}_{i}y_{i,s-l}.$*

Remark 4.2**.**

For model (1), the result in Lemma 4.1 can be applied by having $R=0,L=0,s=1$ and $T=1$. Then, $D_{n+1,s}$ can be written as $D_{n+1}$ containing only $\{p_{1,1},\ldots,p_{n,1}\}$, and we have the conditional distribution

[TABLE]

where $m(D_{n+1})=\alpha_{0}+\mathbf{x}^{\prime}_{n+1}\boldsymbol{\alpha}+\boldsymbol{r}_{\boldsymbol{\theta}}^{\prime}\boldsymbol{R}_{\boldsymbol{\theta}}^{-1}(\log\frac{\boldsymbol{p}_{1}}{\mathbf{1}-\boldsymbol{p}_{1}}-\boldsymbol{\mu}^{n})$, $v(D_{n+1})=\sigma^{2}(1-\boldsymbol{r}_{\boldsymbol{\theta}}^{\prime}\boldsymbol{R}_{\boldsymbol{\theta}}^{-1}\boldsymbol{r}_{\boldsymbol{\theta}}),$ $\boldsymbol{\mu}^{n}=(\alpha_{0}+\mathbf{x}^{\prime}_{1}\boldsymbol{\alpha},\ldots,\alpha_{0}+\mathbf{x}^{\prime}_{n}\boldsymbol{\alpha})^{\prime},\boldsymbol{r}_{\boldsymbol{\theta}}=(R_{\boldsymbol{\theta}}(\mathbf{x}_{n+1},\mathbf{x}_{1}),\ldots,R_{\boldsymbol{\theta}}(\mathbf{x}_{n+1},\mathbf{x}_{n}))^{\prime}$, and $\boldsymbol{R}_{\boldsymbol{\theta}}=\{R_{\boldsymbol{\theta}}(\mathbf{x}_{i},\mathbf{x}_{j})\}$.

Based on Lemma 4.1, the prediction of $p_{s}(\mathbf{x}_{n+1})$ for some untried setting $\mathbf{x}_{n+1}$ and its variance can then be obtained in the next theorem. The proof is given in the supplementary material S6. The definition of minimum mean squared prediction error of $p$ given $D$ is first stated as follows,

[TABLE]

where $F(\cdot)$ is the joint distribution of $(p,D)$ and $\mathbb{E}_{F}[\cdot]$ denotes expectation under distribution $F(\cdot)$.

Theorem 4.3**.**

Given $D_{n+1,s}=\{\boldsymbol{y}^{\prime}_{1},\ldots,\boldsymbol{y}^{\prime}_{T},\boldsymbol{p}^{\prime}_{1},\ldots,\boldsymbol{p}^{\prime}_{T},y_{n+1,s-1},\ldots,y_{n+1,1},p_{n+1,s-1},\ldots,p_{n+1,1}\}$,

• (i)

the minimum mean squared prediction error (MMSPE) predictor of $p_{s}(\mathbf{x}_{n+1})$, denoted by $\hat{p}_{s}(\mathbf{x}_{n+1})$, is

[TABLE]

• (ii)

the MMSPE predictor is an interpolator, i.e., if $\mathbf{x}_{n+1}=\mathbf{x}_{i}$ for $i=1,\cdots,n$, then $\hat{p}_{s}(\mathbf{x}_{n+1})=\mathbb{E}\left[p_{s}(\mathbf{x}_{n+1})|D_{n+1,s}\right]=p_{s}(\mathbf{x}_{i})$ and the predictive variance is 0;

• (iii)

the $q$-th quantile of the conditional distribution $p(\mathbf{x}_{n+1})|D_{n+1,s}$ is

[TABLE]

where $z_{q}$ is the $q$-th quantile of the standard normal distribution.

Theorem 4.3 shows that, given $D_{n+1,s}$, the new predictor for binary data can interpolate the underlying probabilities which generate the training data. According to Theorem 4.3(iii) and the fact that $v(D_{n+1,s})$ increases with the distance to the training data, this result shows an increasing predictive uncertainty for points away from the training data. This predictive property is desirable and consistent with the conventional GP predictor.

In practice, only the binary outputs are observable and the underlying probabilities are not available in the training data. Thus, the following results construct the MMSPE predictor of $p_{s}(\mathbf{x}_{n+1})$ given $\boldsymbol{Y}=(\boldsymbol{y}^{\prime}_{1},\ldots,\boldsymbol{y}^{\prime}_{T},y_{1}(\mathbf{x}_{n+1}),\ldots,y_{s-1}(\mathbf{x}_{n+1}))^{\prime}$. These results can be used for prediction and quantification of the predictive uncertainty, such as constructing predictive confidence intervals for untried settings.

Theorem 4.4**.**

Given $\boldsymbol{Y}=(\boldsymbol{y}^{\prime}_{1},\ldots,\boldsymbol{y}^{\prime}_{T},y_{1}(\mathbf{x}_{n+1}),\ldots,y_{s-1}(\mathbf{x}_{n+1}))^{\prime}$,

• (i)

The MMSPE predictor of $p_{s}(\mathbf{x}_{n+1})$ is

[TABLE]

with variance

[TABLE]

where $\boldsymbol{p}=(\boldsymbol{p}^{\prime}_{1},\cdots,\boldsymbol{p}^{\prime}_{T},p_{1}(\mathbf{x}_{n+1}),\ldots,p_{s-1}(\mathbf{x}_{n+1}))^{\prime}$.

• (ii)

When $\mathbf{x}_{n+1}=\mathbf{x}_{i}$, the MMSPE predictor becomes $\hat{p}_{s}(\mathbf{x}_{i})=\mathbb{E}_{\boldsymbol{p}|\boldsymbol{Y}}\left[p_{s}(\mathbf{x}_{i})|\boldsymbol{Y}\right]$ with variance $\mathbb{V}_{\boldsymbol{p}|\boldsymbol{Y}}\left[p_{s}(\mathbf{x}_{i})|\boldsymbol{Y}\right]$.

• (iii)

The quantiles of $p_{s}(\mathbf{x}_{n+1})|\boldsymbol{Y}$ can be approximated by the sample quantiles of $\{p^{(j)}_{s}\}$, where $p^{(j)}_{s}$ are random samples generated from $p_{s}(\mathbf{x}_{n+1})|D_{n+1,s}$ following the logit-normal distribution given in Lemma 4.1, and $\boldsymbol{p}$ in $D_{n+1,s}=\{\boldsymbol{p},\boldsymbol{Y}\}$ are random samples from distribution $\boldsymbol{p}|\boldsymbol{Y}$.

Although the outcomes in the study are binary, the major focus is on the predictive uncertainty of the underlying probability $p_{s}(\mathbf{x}_{n+1})$. There are two reasons. First, it is more informative and usually of scientific interest to predict generating probabilities instead of binary outcomes. Second, given the understanding of predictive distribution for $p_{s}(\mathbf{x}_{n+1})$ in Theorem 4.4, the predictive uncertainty of $y_{s}(\mathbf{x}_{n+1})$ can be easily constructed by the use of bootstrap predictive distribution, which is given in the following remark.

Remark 4.5**.**

The bootstrap predictive distribution of $y_{s}(\mathbf{x}_{n+1})$ follows $Bernoulli(p_{s}(\mathbf{x}_{n+1}))$, where the distribution of $p_{s}(\mathbf{x}_{n+1})$ can be obtained from Theorem 4.4 (iii). The detailed algorithm can be found in the supplementary material S7.

Remark 4.6**.**

For model (1), the results of Theorem 4.3 and Theorem 4.4 can be applied by assuming $s=1$ and $T=1$.

Although there is no closed form expression for the distribution of $\boldsymbol{p}|\boldsymbol{Y}$, the random samples from $\boldsymbol{p}|\boldsymbol{Y}$ can be easily generated by the Metropolis-Hastings (MH) algorithm. For example, let $\{\boldsymbol{p}^{(j)}\}_{j=1,\ldots,J}$ be the $J$ random samples generated from distribution $\boldsymbol{p}|\boldsymbol{Y}$, then the MMSPE predictor of $p_{s}(\mathbf{x}_{n+1})$ in Theorem 4.4 can be approximated using Monte Carlo method by

[TABLE]

where $D^{(j)}_{n+1,s}=\{\boldsymbol{p}^{(j)},\boldsymbol{Y}\}$. Similar idea can be applied to compute $\mathbb{V}\left[p_{s}(\mathbf{x}_{n+1})|\boldsymbol{Y}\right]$ and the predictive quantiles. Details are given in the supplementary material S7.

Without the information of the underlying probabilities, the predictor does not interpolate all the training data as in Theorem 4.3 (ii). From Theorem 4.4, when $\mathbf{x}_{n+1}=\mathbf{x}_{i}$, the predictor is still unbiased but the corresponding variance is nonzero. Instead, the variance becomes $\mathbb{V}_{\boldsymbol{p}|\boldsymbol{Y}}\left[p_{s}(\mathbf{x}_{i})|\boldsymbol{Y}\right]$, which is due to the uncertainty of the underlying probability. The proof is similar to Theorem 4.3 (ii). To show the empirical performance of the predictive distribution in Theorem 4.4, a one-dimensional example is illustrated in Figure 1. Consider the true probability function, $p(x)=0.4\exp(-1.2x)\cos(3.5\pi x)+0.4$, which is represented by a black dotted line, and the training set that contains 12 evenly-spaced inputs and the corresponding binary outputs represented by red dots. The blue line is the MMSPE predictor constructed by equation (13) and the gray region is the corresponding 95% confidence band constructed by the 2.5%- and 97.5%-quantiles. It appears that the proposed predictor and the confidence band reasonably capture the underlying probability.

When the historical time series for an untried setting (i.e., $y_{1}(\mathbf{x}_{n+1}),\ldots,y_{s-1}(\mathbf{x}_{n+1})$ in Theorem 4.4) is not available, we can emulate a completely new time series (or batch of time series) with input $\mathbf{x}_{n+1}$. The idea is to generate draws from the conditional distribution $p_{s}(\mathbf{x}_{n+1})|\boldsymbol{Y}$ for future outputs, starting from $s=0$, and take pointwise median of the random draws. This idea is similar to the dynamic emulators introduced by Liu and West, (2009) for continuous outputs. The random samples from $p_{s}(\mathbf{x}_{n+1})|\boldsymbol{Y}$ can be generated by the fact $f(p_{s}(\mathbf{x}_{n+1}),\boldsymbol{p}|\boldsymbol{Y})=f(\boldsymbol{p}|\boldsymbol{Y})f(p_{s}(\mathbf{x}_{n+1})|\boldsymbol{p},\boldsymbol{Y})$, where $f(p_{s}(\mathbf{x}_{n+1})|\boldsymbol{p},\boldsymbol{Y})$ is a logit-normal distribution provided in Lemma 4.1. As mentioned above, the random samples from $f(\boldsymbol{p}|\boldsymbol{Y})$ can be generated through the MH algorithm. Therefore, generating a draw from $p_{s}(\mathbf{x}_{n+1})|\boldsymbol{Y}$ consists of two steps: (i) generating the “previous” probability values $\boldsymbol{p}^{*}$ given output $\boldsymbol{Y}$ from the distribution $\boldsymbol{p}|\boldsymbol{Y}$ through the MH algorithm, and (ii) based on the sample $\boldsymbol{p}^{*}$, draw a sample $p^{*}_{s}(\mathbf{x}_{n+1})$ from $p_{s}(\mathbf{x}_{n+1})|\boldsymbol{p}^{*},\boldsymbol{Y}$, which is a logit-normal distribution, and also draw a sample $y^{*}_{s}(\mathbf{x}_{n+1})$ from a Bernoulli distribution with parameter $p^{*}_{s}(\mathbf{x}_{n+1})$. An explicit algorithm is given in the supplementary material S8.

5 Simulation Studies

In Section 5.1, we conduct simulations generated from Gaussian processes to demonstrate the estimation performance. In Section 5.2, the prediction performance is examined by comparing several existing methods using the data generated from a modified Friedman function (Friedman,, 1991).

5.1 Estimation Performance

Consider a 5-dimensional input space, $d=5$, and the input $\mathbf{x}$ is randomly generated from a regular grid on $[0,1]^{5}$. The binary output, $y_{t}(\mathbf{x})$ at time $t$, is simulated by a Bernoulli distribution with probability $p_{t}(\mathbf{x})$ calculated by (3) and $\alpha_{0}=0.5$, $\boldsymbol{\alpha}=(-3,2,-2,1,0.5)^{\prime}$, $\varphi_{1}=0.8$, $\sigma^{2}=1$, and the power exponential correlation function (2) is chosen with $\boldsymbol{\theta}=(0.5,1.0,1.5,2.0,2.5)^{\prime}$ and $p=2$. Four sample size combinations of $n$ and $T$ are considered in the simulations.

The potential confounding between the polynomials in the mean function and the zero-mean Gaussian process can lead to the lack of identifiability, which will cause the estimated mean model to lose interpretability. In order to tackle this problem, Plumlee and Joseph, (2018) proposed an orthogonal Gaussian process model whose stochastic part is orthogonal to the mean function. The key idea is to construct the correlation function that achieves the orthogonality. The orthogonal correlation function derived from the exponential correlation functions with power $p=2$ is given in equation (8) of Plumlee and Joseph, (2018). We implemented the orthogonal correlation function (abbreviated as OGP) as well as the power exponential correlation function (abbreviated as PE) in the simulation.

The estimation results for the linear function coefficients are summarized in Table 1 based on 100 replicates for each sample size combination. In general, the proposed approach can estimate the linear function coefficients ($\alpha_{0},\boldsymbol{\alpha},\varphi_{1}$) reasonably well. Compared with PE, the estimation improvement using OGP is reported as IMP. It appears that the estimation accuracy can be further improved by the use of orthogonal correlation functions. Therefore, orthogonal correlation functions are generally recommended when estimation is of major interest, such as in variable selection and calibration problems.

The parameter estimation results for $\sigma^{2}$ and PE correlation parameters $\boldsymbol{\theta}$ are reported in Table 2. The estimation with OGP has similar results, so we omit them to save space. The proposed approach tends to overestimate the correlation parameters for small sample size. This is not surprising because the estimation of correlation parameters is more challenging and the same phenomenon is observed in conventional GP models (see Li and Sudjianto, (2005)). This problem can be ameliorated by the increase of sample size as shown in Table 2. Given the same number of total sample size, ($n=200,T=50$) and ($n=500,T=20$), it appears that a larger $n$ can improve the estimation accuracy more effectively.

Based on the construction of predictive distribution in Section 4, we can emulate a new time series with an untried input. Here we generate 100 random untried inputs to examine its prediction performance. The prediction performance is evaluated by the following two measures. Define the 100 random untried inputs ($n_{\text{test}}=100$) by $\mathbf{x}^{*}_{1},\ldots,\mathbf{x}^{*}_{100}$. Since the underlying probabilities are known in the simulation settings, we evaluate the prediction performance by the root mean squared prediction error

[TABLE]

where $\hat{p}_{t}(\mathbf{x}^{*}_{i})$ is the predictive probability. The results are given in Table 3. Overall, the proposed predictor has the root mean squared prediction error less than $0.12$. Also, with the increase of sample size, the prediction error decreases in general.

Furthermore, the predictive distributions can be used to quantify the prediction uncertainty. The predictive distributions with two random untried inputs are shown in Figure 2, where the green dotted lines represent the true probability, the red dashed lines represent the MMSPE predictors obtained in Theorem 4.4. From Figure 2, it appears that the MMSPE predictors provide accurate predictions in both cases. Moreover, the predictive distributions provide rich information for statistical inference. For example, we can construct $95\%$ predictive confidence intervals for the two untried settings as indicated in blue in Figure 2.

5.2 Prediction Performance

To examine the performance of the proposed model as an emulator, we compare its prediction accuracy with four existing methods: (1) the logistic regression model, (2) a combination of logistic regression model with time series mean function, (3) the Bayesian generalized Gaussian process model (Williams and Barber,, 1998), which incorporates a Gaussian process prior but does not take into account the time series structure, and (4) the functional Gaussian process proposed by Shi and Choi, (2011), which captures the serial correlation by functional data analysis techniques. These methods are respectively implemented by R (R Core Team,, 2015) using packages binaryGP (Sung,, 2017), stat (R Core Team,, 2015), a modification of stat, kernlab (Karatzoglou et al.,, 2004) and GPFDA (Shi and Cheng,, 2014) adapted to classification.

The simulated data are generated by a modification of the Friedman function (Friedman,, 1991),

[TABLE]

where $\mathbf{x}\in[0,1]^{5}$ and the Friedman function is given in the brackets with intercept $-5$ and scale $1/3$ to ensure $p_{t}(\mathbf{x})$ is uniformly located at $[0,1]$. The input $\mathbf{x}$ is randomly generated from $[0,1]^{5}$ and the corresponding binary output $y_{t}(\mathbf{x})$ at time $t$ is generated by a Bernoulli distribution with probability $p_{t}(\mathbf{x})$. The size of the training data is set to be $n=200,T=20$.

Since the underlying probabilities are known in this simulation setting, the prediction performance is evaluated by RMSPE using 100 randomly generated untried settings $(n_{\rm test}=100,T=20)$. The results for the five methods based on 100 replicates are shown in the left panel of Figure 3. In general, the proposed method has lower RMSPE than the other four methods. By incorporating a Gaussian process to model the nonlinearity, the proposed method outperforms the straightforward combination of logistic regression model and time series structure. On the other hand, comparing with Bayesian generalized Gaussian process model (i.e., kernlab in Figure 3), the proposed method further improves the prediction accuracy by taking into account the time series structure. The computation time, for model fitting and prediction, is given in the right panel of Figure 3. The proposed method is faster than GPFDA. Comparing with glm and glm_ts, the major computational difficulty lies in the estimation of correlation parameters, which is a common issue in conventional GPs because there is no analytical solution for the parameter estimation. The kernlab has better computational performance since it assumes that all the correlation parameters are the equal and estimated by analytic approximation, that is, $\theta_{1}=\ldots,\theta_{d}$ in (2). However, the computation time of kernlab is expected to increase if this assumption is relaxed to a correlation function with different correlation parameters as in (2).

6 Computer Experiments for Cell Adhesion Frequency Assay

In an earlier study based on in vitro experiments, an important memory effect was discovered in the repeated adhesion experiments of the micropipette adhesion frequency assay. However, only limited variables of interest can be studied in the lab because of the technical complexity of the biological setting and the complicated experimental manipulation. Therefore, computer simulation experiments are performed to examine the complex mechanisms behind repeated receptor-ligand binding to rigorously elucidate the biological mechanisms behind the memory effect.

In these computer experiments, two surfaces are simulated to reflect the two opposing membranes in the adhesion frequency assays. The molecules on the surfaces are permitted to interact for the contact duration and then separated for a period of waiting time to simulate the retract-approach phase of the assays. The computer experiments are constructed based on a kinetic proofreading model for receptor modification and solved through a Gillespie algorithm (Gillespie,, 1976), which is a stochastic simulation algorithm. The contact is scored as 1 or 0 depending on whether at least one bond or no bond is observed, respectively. The process is repeated until the given number of contacts is completed.

The biological system investigated here is the T Cell Receptor (TCR) binding to antigen peptide bound to Major Histocompatibility Complex (pMHC), which has previously been shown to exhibit memory in repeated contacts (Zarnitsyna et al.,, 2007). The TCR is the primary molecule involved in detecting foreign antigens which are presented on pMHC molecules expressed by infected cells. Memory in serial interactions of these foreign antigens may be a mechanism which underlies the major properties of T cell antigen recognition: sensitivity, specificity, and context discrimination. It has largely remained uninvestigated due to the small time scales at which the mechanism operates and the complexity of the experimental system. Although there are many possible cellular mechanisms which may induce this behavior, we investigate a specific mechanism, called free induction mechanism (Huang et al.,, 2010), in this study as to how this memory may be controlled: pMHC binding to a single TCR within a cluster upregulates the kinetics of all TCRs within that cluster.

The free induction mechanism has six control variables given in Table 4. The range of each control variable in Table 4 is given by physical principles or estimated through similar molecular interactions. The design of the control variables is a 60-run OA-based Latin hypercube designs (Tang,, 1993). For each run, it consists of 50 replicates and each replicate has 100 repeated contacts ($T=100$).

The proposed estimation method is implemented with orthogonal correlation functions derived from the power exponential correlation function with $p=2$, which can be found in equation (8) in Plumlee and Joseph, (2018). We start with a large model in which the mean function includes all the main effects of the control variables and their interactions with the past time series output $y_{t-1}$. The model is written as:

[TABLE]

where all the control variables are standardized to $[0,1]$, $\hat{\sigma}=0.43$ and the estimated correlation parameters are $\hat{\boldsymbol{\theta}}=(\hat{\theta}_{K_{f,p}},\hat{\theta}_{K_{r,p}},\hat{\theta}_{T_{half}},\hat{\theta}_{Tc},\hat{\theta}_{Tw},\hat{\theta}_{Kc})=(3.28,1.70,7.77,0.06,4.78,0.74)$. Estimation results for the mean function coefficients are given in Table 5 with p values calculated based on the asymptotic results in Theorem 3.1. We use these p values to perform variable selection and identify significant effects for the mean function. According to Table 5, ${x}_{T_{half}}$ has no significant effect in the mean function at the $0.01$ level. By removing ${x}_{T_{half}}$, the model can be updated as

[TABLE]

where $\hat{\sigma}=0.44$, the estimated correlation parameters are $\hat{\boldsymbol{\theta}}=(\hat{\theta}_{K_{f,p}},\hat{\theta}_{K_{r,p}},\hat{\theta}_{T_{half}},\\ \hat{\theta}_{Tc},\hat{\theta}_{Tw},\hat{\theta}_{Kc})=(3.27,1.71,7.77,0.06,4.81,0.74)$. Based on the updated model, among the interaction effects, all the control variables except ${x}_{T_{half}}$ are significant for inducing memory in the free induction model.

The application of this statistical approach to the analysis of simulations and experimental data will be powerful in illuminating the unknown biological mechanism, and also informs the next round of experiments by advising future manipulations. Additionally, developments on the calibration of computer experiments based upon the proposed predictive distribution will help provide insight into the range of possible values of variables, such as the increases in kinetic rates, which are difficult to determine through existing methods due to the small time scale at which this mechanism operates and the limits of existing experimental techniques.

Besides estimation, the proposed method also provides predictors which can serve as efficient and accurate emulators for untried computer experiments. The construction of emulator is an important step for future research on calibration where computer experiment outputs under the same settings of the lab experiments are required but not necessarily available. To assess the predictive performance, we compare the proposed method with the four existing methods discussed in Section 5 based on a 10-fold cross-validation study. To evaluate the binary predictive performance, we consider four proper scoring rules defined in Table 1 of Gneiting and Raftery, (2007), including the Brier score, the spherical score, the logarithmic score, and the zero-one score. The results are reported in Table 6, where larger values indicate better predictions. Compared to the four existing methods, the proposed method consistently has better prediction performance across the four scoring rules. The GPFDA approach appears to be the second best, but it is usually time-consuming to evaluate as we observed in Figure 3. In addition, the 10-fold variation is relatively large as shown in Figure 4 using the zero-one scoring rule.

7 Summary and Concluding Remarks

In spite of the prevalence of Gaussian process models in the analysis of computer experiments, their applications are limited to the Gaussian assumption on the responses. Motivated by the study of cell adhesion where the computer simulation responses are binary time series, a generalized Gaussian process model is proposed in this paper. The estimation procedure is introduced and asymptotic properties are derived. An optimal predictor and its predictive distribution are constructed which can be used for uncertainty quantification and calibration of future computer simulations. An R package is available for implementing the proposed methodology. The methodology is applied to analyze stochastic computer simulations for a cell adhesion mechanism. The results reveal important biological information which is not available in lab experiments and provide valuable insights on how the next round of lab experiments should be conducted.

The current work can be extended in several directions. First, we will extend the proposed method to other non-Gaussian data, such as the count data. It is conceivable that the current estimation procedure can be directly extended to other exponential family distributions, but different predictive distributions are expected for different types of non-Gaussian data. Second, the computational cost in the proposed procedure can be further reduced. In particular, the inversion of $\boldsymbol{R}_{\boldsymbol{\theta}}$ can be computationally prohibitive when sample size is large. This computational issue has been addressed for conventional GP models in the recent literature. Extensions of these methods (e.g., Gramacy and Apley, (2015); Sung et al., (2018)) to binary responses deserve further attention. Third, many mathematical models underlying the computer simulations contain unknown parameters, which need to be estimated using data from lab experiments. This problem is called calibration and much work has been done in the computer experiment literature. However, the existing methods (e.g., Kennedy and O’Hagan, (2001), Tuo and Wu, (2015) and Gramacy et al., (2015)) are only applicable under the Gaussian assumption. Based upon the model and prediction procedure proposed herein, we will work on developing a calibration method for non-Gaussian data.

Supplementary Materials The assumptions for Theorems 3.1 and 3.3, the proofs of Theorems 3.1, 3.3, 4.1, 4.3, and the algorithms for estimation and emulation are given in an online supplement.

Acknowledgements: The authors gratefully acknowledge helpful advice from the associate editor and two referees. This work was supported by NSF DMS 1660504 and 1660477.

Appendix A Algorithm: Estimation of $(\boldsymbol{\beta},\boldsymbol{\omega}$)

1:Set initial values $\boldsymbol{\omega}=(\sigma^{2},\boldsymbol{\theta})=\mathbf{1}_{d+1},\boldsymbol{\beta}=\mathbf{1}_{m},p_{it}=1$, and set $\tilde{\eta}_{it}=\log\frac{p_{it}}{1-p_{it}}+\frac{y_{it}-p_{it}}{p_{it}(1-p_{it})}$ for each $i$ and $t$.

2:repeat

3:     repeat

4:         Set $\boldsymbol{W}$ as an $N\times N$ diagonal matrix with diagonal elements $W_{it}=p_{it}(1-p_{it})$

5:         Set $\boldsymbol{V}=\boldsymbol{W}^{-1}+\sigma^{2}(\boldsymbol{R}_{\boldsymbol{\theta}}\otimes I_{T})$

6:         Update $\boldsymbol{\beta}=(\boldsymbol{X}^{\prime}\boldsymbol{V}^{-1}\boldsymbol{X})^{-1}\boldsymbol{X}^{\prime}\boldsymbol{V}^{-1}\tilde{\boldsymbol{\eta}}$

7:         Set $\boldsymbol{Z}=\sigma^{2}(\boldsymbol{R}_{\boldsymbol{\theta}}\otimes I_{T})\boldsymbol{V}^{-1}(\tilde{\boldsymbol{\eta}}-\boldsymbol{X}^{\prime}\boldsymbol{\beta})$

8:         Update $p_{it}=\left(\frac{\exp\{\boldsymbol{X}^{\prime}\boldsymbol{\beta}+\boldsymbol{Z}\}}{\mathbf{1}_{N}+\exp\{\boldsymbol{X}^{\prime}\boldsymbol{\beta}+\boldsymbol{Z}\}}\right)_{it}$ and $\tilde{\eta}_{it}=\log\frac{p_{it}}{1-p_{it}}+\frac{y_{it}-p_{it}}{p_{it}(1-p_{it})}$ for each $i$ and $t$

9:     until $\{\tilde{\eta}_{it}\}_{it}$ converges

10:     Update $\boldsymbol{\omega}=\arg\min_{\boldsymbol{\omega}}L(\boldsymbol{\omega})$, where $L(\boldsymbol{\omega})$ is the negative log-likelihood function (12)

11:     Update $(\sigma^{2},\boldsymbol{\theta})=\boldsymbol{\omega}$

12:until $\boldsymbol{\beta}$ and $\boldsymbol{\omega}$ converge

13:Return $\boldsymbol{\beta}$ and $\boldsymbol{\omega}$

Appendix B Assumptions

The parameter $\boldsymbol{\beta}$ belongs to an open set $B\subseteq\mathbb{R}^{m}$ and the parameter $\boldsymbol{\omega}$ belongs to an open set $\Omega\subseteq\mathbb{R}^{d+1}$. 2. 2.

The model matrix $X_{it}$ lies almost surely in a nonrandom compact subset of $\mathbb{R}^{m}$ such that $Pr(\sum^{n}_{i=1}\sum^{T}_{t=1}X^{\prime}_{it}X_{it}>0)=1$.

For any matrix $A$, define $\|A\|\equiv\sqrt{\text{tr}(A^{\prime}A)}$; for the covariance matrix $\boldsymbol{V}(\boldsymbol{\omega})$, define $V_{i}(\boldsymbol{\omega})\equiv\partial\boldsymbol{V}(\boldsymbol{\omega})/\partial\omega_{i}$ and $V_{ij}(\boldsymbol{\omega})\equiv\partial\boldsymbol{V}(\boldsymbol{\omega})/\partial\omega_{i}\partial\omega_{j}$; for $\boldsymbol{\omega}\in\Omega$, denote $\stackrel{{\scriptstyle u}}{{\longrightarrow}}$ as uniform convergence of nonrandom functions over compact subsets of $\Omega$.

$J_{N}(\boldsymbol{\omega})P_{N}(\boldsymbol{\omega})^{-1}\stackrel{{\scriptstyle d}}{{\longrightarrow}}W(\boldsymbol{\omega})$ for some nonsingular $W(\boldsymbol{\omega})$, which is continuous in $\boldsymbol{\omega}$, where $P_{N}(\boldsymbol{\omega})=\text{diag}(\|\varPi(\boldsymbol{\omega})V_{1}(\boldsymbol{\omega})\|,\ldots,\|\varPi(\boldsymbol{\omega})V_{d+1}(\boldsymbol{\omega})\|)$ and $\varPi(\boldsymbol{\omega})=\boldsymbol{V}(\boldsymbol{\omega})^{-1}-\boldsymbol{V}(\boldsymbol{\omega})^{-1}\boldsymbol{X}(\boldsymbol{X}^{\prime}\boldsymbol{V}(\boldsymbol{\omega})^{-1}\boldsymbol{X})^{-1}\boldsymbol{X}^{\prime}\boldsymbol{V}(\boldsymbol{\omega})^{-1}$. 2. 4.

If there exists a sequence $\{r_{N}\}_{N\geq 1}$ with $\lim\sup_{N\rightarrow\infty}r_{N}/N\leq 1-\delta$, for some $\delta\in(0,1)$, such that for any compact subset $K\subseteq\Omega$, there exist constants $0<C_{1}(K)<\infty$ and $C_{2}(K)>0$ such that

[TABLE]

and

[TABLE]

uniformly in $\boldsymbol{\omega}\in K$, where $|\lambda_{1}|\leq\ldots\leq|\lambda_{N}|$ are the absolute eigenvalues of $\boldsymbol{V}(\boldsymbol{\omega})$, $|\lambda^{i}_{1}|\leq\ldots\leq|\lambda^{i}_{N}|$ are the absolute eigenvalues of $V_{i}(\boldsymbol{\omega})$, and $|\lambda^{ij}_{1}|\leq\ldots\leq|\lambda^{ij}_{N}|$ are the absolute eigenvalues of $V_{ij}(\boldsymbol{\omega})$.

Assumption 2 holds when the row vectors of $\boldsymbol{X}$ are linear independent. Thus, if only the linear effect is considered in the mean function, then orthogonal designs or orthogonal array-based designs, such as OA-based Latin hypercube designs (Tang,, 1993), can be chosen for sampling schemes. The conditions for Assumption 4 can be referred to Cressie and Lahiri, (1996), in which the checkable conditions for rectangular lattice of data sites and irregularly located data sites are given. For instance, for rectangular lattice of data sites, with certain correlation functions, a sufficient condition is choosing data locations whose minimum distance is sufficiently large. More details can be seen in Cressie and Lahiri, (1996). Thus, space-filling designs, such as Latin hypercube designs (McKay et al.,, 1979) and maximin distance designs (Johnson et al.,, 1990), can be chosen for sampling schemes.

Appendix C Proof of Theorem 3.1

The model (4) can be seen as a binary time series model with random effects by multiplying an identity matrix on $\boldsymbol{Z}$, that is,

[TABLE]

where $I_{N}$ and $\boldsymbol{Z}$ are viewed as the model matrix and coefficients of random effects, respectively. Therefore, if the variance-covariance parameters are given, inference of $\boldsymbol{\beta}$ is a special case of the binary time series model with random effects in Hung et al., (2008). Therefore, following Theorem 1 in Hung et al., (2008), the score function $S_{N}(\boldsymbol{\beta},\boldsymbol{\omega})$ is asymptotically normally distributed.

Appendix D Proof of Theorem 3.3

According to Breslow and Clayton, (1993), one can view the inference on the variance-variance component as an iterative procedure for the linear mixed model

[TABLE]

with the iterative weight $\boldsymbol{W}^{-1}$. Thus, it is a special case of the Gaussian general linear model in Cressie and Lahiri, (1993) with response vector $\tilde{\boldsymbol{\eta}}$ and variance-covariance component $\Sigma(\boldsymbol{\omega})+\boldsymbol{W}^{-1}$ with parameters $\boldsymbol{\omega}$. Since the asymptotic distribution of REML estimators for the variance-covariance parameters has been shown in Cressie and Lahiri, (1993) for a Gaussian general linear model, the result directly follows as a special case of Corollary 3.3 in Cressie and Lahiri, (1993). Note that Assumption 4 in the supplementary material B implies the conditions for Corollary 3.3 in Cressie and Lahiri, (1993). See the proof of Theorem 2.2 in Cressie and Lahiri, (1996).

Appendix E Proof of Lemma 4.1

We start the proof by deriving the conditional distribution from a simple model (1) (without time-series), and then extend the result to prove Lemma 4.1. First, a definition and a lemma about multivariate log-normal distribution are in order.

Definition E.1**.**

Suppose $\boldsymbol{\xi}=(\xi_{1},\ldots,\xi_{n})^{\prime}$ has a multivariate normal distribution with mean $\boldsymbol{\mu}_{n}$ and covariance variance $\boldsymbol{\Sigma}_{n\times n}$. Then $\boldsymbol{b}=\exp\{\boldsymbol{\xi}\}$ has a multivariate log-normal distribution. Denote it as $\boldsymbol{b}\sim\mathcal{LN}(\boldsymbol{\mu}_{n},\boldsymbol{\Sigma}_{n\times n})$.

Lemma E.1**.**

Suppose $\boldsymbol{b}^{n}$ and $b_{n+1}$ have a multivariate log-normal distribution

[TABLE]

The conditional distribution of $b_{n+1}$ given $\boldsymbol{b}^{n}$ is $b_{n+1}|\boldsymbol{b}^{n}\sim\mathcal{LN}(\mu^{*},v^{*})$, where $\mu^{*}=\mu_{n+1}+\boldsymbol{r}^{\prime}\boldsymbol{\Sigma}_{n\times n}^{-1}(\log\boldsymbol{b}^{n}-\boldsymbol{\mu}^{n})$ and $v^{*}=\sigma^{2}_{n+1}-\boldsymbol{r}^{\prime}\boldsymbol{\Sigma}_{n\times n}^{-1}\boldsymbol{r}.$

Proof.

Using transformation of a standard normal distribution, one can show that the joint probability density function of the multivariate log-normal distribution $\boldsymbol{b}^{n}$ is

[TABLE]

Denote $\boldsymbol{b}^{n+1}=(b_{1},\ldots,b_{n},b_{n+1})$, $\boldsymbol{\mu}^{n+1}=(\mu_{1},\ldots,\mu_{n},\mu_{n+1})$ and

[TABLE]

Then, the conditional probability density function of $b_{n+1}$ given $\boldsymbol{b}^{n}$ can be derived as

[TABLE]

Let $\boldsymbol{a}_{1}=\log\boldsymbol{b}^{n}-\boldsymbol{\mu}^{n}$ and $\boldsymbol{a}_{2}=\log b^{n+1}-\mu^{n+1}$. Applying the partitioned matrix inverse results (page 99 of Harville, (1997)) gives

[TABLE]

where $\sigma_{22\cdot 1}=\sigma^{2}_{n+1}-\boldsymbol{r}^{\prime}\boldsymbol{\Sigma}^{-1}_{n\times n}\boldsymbol{r}$ and is a real number.

Thus, the conditional probability density function of $b_{n+1}$ given $\boldsymbol{b}^{n}$ can be simplified as

[TABLE]

Therefore, according to the probability density function of a log-normal distribution, we have $b_{n+1}|\boldsymbol{b}^{n}\sim\mathcal{LN}(\mu^{*},v^{*})$, where $\mu^{*}=\mu_{n+1}+\boldsymbol{r}^{\prime}\boldsymbol{\Sigma}^{-1}_{n\times n}(\log\boldsymbol{b}^{n}-\boldsymbol{\mu}^{n})$ and $v^{*}=\sigma_{22\cdot 1}=\sigma^{2}_{n+1}-\boldsymbol{r}^{\prime}\boldsymbol{\Sigma}_{n\times n}^{-1}\boldsymbol{r}$. ∎

Lemma E.2**.**

Consider the model (1) (without time-series), given $(p(\mathbf{x}_{1}),\ldots,p(\mathbf{x}_{n}))^{\prime}=\boldsymbol{p}^{n}$, the conditional distribution of $p(\mathbf{x}_{n+1})$ is a logit-normal distribution, that is, $p(\mathbf{x}_{n+1})|\boldsymbol{p}^{n}\sim Logitnormal(m(\boldsymbol{p}^{n}),v(\boldsymbol{p}^{n}))$ with

[TABLE]

where $\boldsymbol{\mu}^{n}=(\mu(\mathbf{x}_{1}),\ldots,\mu(\mathbf{x}_{n}))^{\prime}$, $\mu(\mathbf{x}_{i})=\alpha_{0}+\mathbf{x}^{\prime}_{i}\boldsymbol{\alpha}$, $\boldsymbol{r}_{\boldsymbol{\theta}}=(R_{\boldsymbol{\theta}}(\mathbf{x}_{n+1},\mathbf{x}_{1}),\ldots,R_{\boldsymbol{\theta}}(\mathbf{x}_{n+1},\mathbf{x}_{n}))^{\prime}$, and $\boldsymbol{R}_{\boldsymbol{\theta}}=\{R_{\boldsymbol{\theta}}(\mathbf{x}_{i},\mathbf{x}_{j})\}$.

Proof.

Let $\eta_{i}=\mu(\mathbf{x}_{i})+Z(\mathbf{x}_{i})$ and $b_{i}=\exp\{\eta_{i}\}=p(\mathbf{x}_{i})/(1-p(\mathbf{x}_{i}))$ for $i=1,\ldots,n+1$. Since $(\eta_{1},\ldots,\eta_{n},\eta_{n+1})^{\prime}\sim\mathcal{N}(\boldsymbol{\mu}^{n+1},\sigma^{2}\boldsymbol{R}^{*}_{\boldsymbol{\theta}})$, where $\boldsymbol{\mu}^{n+1}=((\boldsymbol{\mu}^{n})^{\prime},\mu(\mathbf{x}_{n+1}))^{\prime}$ and

[TABLE]

we have $(b_{1},\ldots,b_{n},b_{n+1})^{\prime}\sim\mathcal{LN}(\boldsymbol{\mu}^{n+1},\sigma^{2}\boldsymbol{R}^{*}_{\boldsymbol{\theta}})$ by Definition E.1. Thus, using Jacobian of the transformation and Lemma E.1, we have

[TABLE]

Therefore, according to the probability density function of a logit-normal distribution, we have $p(\mathbf{x}_{n+1})|\boldsymbol{p}^{n}\sim Logitnormal(m(\boldsymbol{p}^{n}),v(\boldsymbol{p}^{n}))$. ∎

Similarly, the result of Lemma E.2 can be extended to the general model (3). Given $\boldsymbol{Y}=(\mathbf{y}^{\prime}_{1},\ldots,\mathbf{y}^{\prime}_{T},y_{n+1,1},\ldots,y_{n+1,s-1})^{\prime}$, at a fixed time-step $s$, $p_{s}(\mathbf{x}_{i})$ can be seen to have the model (1) with mean function $\mu(\mathbf{x}_{i},\boldsymbol{Y})=\sum^{R}_{r=1}\varphi_{r}y_{i,s-r}+\alpha_{0}+\mathbf{x}^{\prime}_{i}\boldsymbol{\alpha}+\sum^{L}_{l=1}\boldsymbol{\gamma}_{l}\mathbf{x}_{i}y_{i,s-l}$. Thus, by Lemma E.2, denote $\boldsymbol{p}_{s}=(p_{s}(\mathbf{x}_{1}),\ldots,p_{s}(\mathbf{x}_{n}))^{\prime}$, we have

[TABLE]

where $m(\boldsymbol{p}_{s},\boldsymbol{Y})=\mu(\mathbf{x}_{n+1},\boldsymbol{Y})+\boldsymbol{r}_{\boldsymbol{\theta}}^{\prime}\boldsymbol{R}_{\boldsymbol{\theta}}^{-1}(\log\frac{\boldsymbol{p}_{s}}{\mathbf{1}-\boldsymbol{p}_{s}}-\boldsymbol{\mu}^{n}),\boldsymbol{\mu}^{n}=(\mu(\mathbf{x}_{1},\boldsymbol{Y}),\ldots,\mu(\mathbf{x}_{n},\boldsymbol{Y}))^{\prime}$, and $v(\boldsymbol{p}_{s},\boldsymbol{Y})=\sigma^{2}(1-\boldsymbol{r}_{\boldsymbol{\theta}}^{\prime}\boldsymbol{R}_{\boldsymbol{\theta}}^{-1}\boldsymbol{r}_{\boldsymbol{\theta}}).$ By the fact that $Z_{t}(\mathbf{x})$ is independent over time, which implies $p_{s}(\mathbf{x})$ is independent of $p_{t}(\mathbf{x})$ for any $t\neq s$, $p_{s}(\mathbf{x}_{n+1})|D_{n+1,s}$ and $p_{s}(\mathbf{x}_{n+1})|\boldsymbol{p}_{s},\boldsymbol{Y}$ have the same distribution. So, $p_{s}(\mathbf{x}_{n+1})|D_{n+1,s}\sim Logitnormal(m(D_{n+1,s}),v(D_{n+1,s}))$, where $m(D_{n+1,s})=m(\boldsymbol{p}_{s},\boldsymbol{Y})$ and $v(D_{n+1,s})=v(\boldsymbol{p}_{s},\boldsymbol{Y})$.

Appendix F Proof of Theorem 4.3

(i) First, one can show that if $(p_{s}(\mathbf{x}_{n+1}),D_{n+1,s})$ has a joint distribution for which the conditional mean of $p_{s}(\mathbf{x}_{n+1})$ given $D_{n+1,s}$ exists, then $\mathbb{E}\left[p(\mathbf{x}_{n+1})|D_{n+1,s}\right]$ is the minimum mean squared error predictor of $p(\mathbf{x}_{n+1})$. See Theorem 3.2.1 in Santner et al., (2003). Thus, by the result of Lemma 4.1, we have the conditional mean $\mathbb{E}\left[p(\mathbf{x}_{n+1})|D_{n+1,s}\right]=\kappa(m(D_{n+1,s}),v(D_{n+1,s}))$ with variance $\mathbb{V}\left[p(\mathbf{x}_{n+1})|D_{n+1,s}\right]=\tau(m(D_{n+1,s}),v(D_{n+1,s}))$.

(ii) If $\mathbf{x}_{n+1}=\mathbf{x}_{i}$ for $i=1,\ldots,n$, then $m(D_{n+1,s})=\log(p_{s}(\mathbf{x}_{i})/(1-p_{s}(\mathbf{x}_{i})))$ and $v(D_{n+1,s})=0$, which implies that

[TABLE]

and $\tau(m(D_{n+1,s}),0)=0$ by using transformation of a normal distribution. Thus, by Theorem 4.3 (i), we have $\mathbb{E}\left[p_{s}(\mathbf{x}_{n+1})|D_{n+1,s}\right]=p_{s}(\mathbf{x}_{i})$ and $\mathbb{V}\left[p_{s}(\mathbf{x}_{n+1})|D_{n+1,s}\right]=0.$

(iii) Let $X\sim\mathcal{N}(m(D_{n+1,s}),v(D_{n+1,s}))$, $P=\exp\{X\}/(1+\exp\{X\})$, which has the distribution $Logitnormal(m(D_{n+1,s}),v(D_{n+1,s}))$, and $Q(q;D_{n+1,s})$ be the $q$-th quantile of $P$. Consider the function $f(x)=\log(x/(1-x))$. The derivative is $f^{\prime}(x)=1/(x(1-x))$. Thus, for $0<x<1$ the derivative is positive and the $f(x)$ function is increasing in $x$. Then,

[TABLE]

[TABLE]

Appendix G Metropolis-Hastings Algorithm and Approximation for Theorem 4.4

The Metropolis-Hastings (MH) algorithm for generating random samples from $\boldsymbol{p}|\boldsymbol{Y}$ is given as follows.

1:for $j=1$ to $J$ do

2:     Set $N_{s}=nT+s-1$.

3:     Start with a zero vector $\boldsymbol{p}$ of size $N_{s}$.

4:     for $k=1$ to $N_{s}$ do

5:         Generate a random value $p^{*}_{k}$ from $Logitnormal(m(\boldsymbol{p}_{-k},\mathbf{y}_{-k}),v(\boldsymbol{p}_{-k},\mathbf{y}_{-k}))$.

6:         Generate an uniform random variable $U\sim Unif(0,1)$.

7:         if $U<\min\{1,\frac{f(y_{k}|p^{*}_{k})}{f(y_{k}|p_{k})}\}$ then

8:              Set $\boldsymbol{p}=(p_{1},\ldots,p^{*}_{k},\ldots,p_{N_{s}})$.               

9:     Set $\boldsymbol{p}^{(j)}=\boldsymbol{p}$

10:Return $\{\boldsymbol{p}^{(j)}\}_{j=1,\ldots,J}$.

In the algorithm, we first sample a value for the $k$-th component $p_{k}$ from the conditional distribution of $p_{k}$ given $p_{j},y_{j},j\neq k$, which is $Logitnormal(m(\boldsymbol{p}_{-k},\mathbf{y}_{-k}),v(\boldsymbol{p}_{-k},\mathbf{y}_{-k}))$, where

[TABLE]

in which $\mu_{t}(\mathbf{x}_{i})=\sum^{R}_{r=1}\varphi_{r}y_{i,t-r}+\mathbf{x}^{\prime}_{i}\boldsymbol{\alpha}+\sum^{L}_{l=1}\boldsymbol{\gamma}_{l}\mathbf{x}_{i}y_{i,t-l}$ and $Q_{kj}$ is the $(k,j)$-element of $\boldsymbol{R}^{-1}_{\boldsymbol{\theta}}$. Similar to Zhang, (2002), we use the single-component MH algorithm, that is, to update only a single component at each iteration. Moreover, the proposed distribution $f(p_{k})$ is used for the single MH algorithm, so that the probability of accepting a new $p^{*}_{k}$ is the minimum of 1 and $\frac{f(p^{*}_{k}|y_{k})f(p_{k})}{f(p_{k}|y_{k})f(p^{*}_{k})}\left(=\frac{f(y_{k}|p^{*}_{k})}{f(y_{k}|p_{k})}\right)$.

Based on the samples $\{\boldsymbol{p}^{(j)}\}_{j=1,\ldots,J}$, the mean, variance, and $q$-quantile of $p_{s}(\mathbf{x}_{n+1})|\boldsymbol{Y}$ can be respectively approximated by

[TABLE]

[TABLE]

and the $q$-quantile of $\{p^{(j)}_{s}\}^{J}_{j=1}$, where $p^{(j)}_{s}$ is generated from $Logitnormal(m(\boldsymbol{p}^{(j)},\boldsymbol{Y}),v(\boldsymbol{p}^{(j)},\boldsymbol{Y}))$. Similarly, the distribution of $y_{s}(\mathbf{x}_{n+1})|\boldsymbol{Y}$ can be approximated by the sample distribution of $\{y^{(j)}_{s}\}^{J}_{j=1}$, where $y^{(j)}_{s}$ is generated from a Bernoulli distribution with probability $p^{(j)}_{s}$.

Appendix H Algorithm: Dynamic Binary Emulator

1:for $j=1$ to $J$ do

2:     Set $N=nT$.

3:     Start with a zero vector $\boldsymbol{p}$ of size $N$.

4:     for $i=1$ to $N$ do

5:         Generate a random value $p^{*}_{k}$ from $Logitnormal(m(\boldsymbol{p}_{-k},\mathbf{y}_{-k}),v(\boldsymbol{p}_{-k},\mathbf{y}_{-k}))$.

6:         Generate an uniform random variable $U\sim Unif(0,1)$.

7:         if $U<\min\{1,\frac{f(y_{k}|p^{*}_{k})}{f(y_{k}|p_{k})}\}$ then

8:              Set $\boldsymbol{p}=(p_{1},\ldots,p^{*}_{k},\ldots,p_{N})$.               

9:     Set $\boldsymbol{p}_{n+1}=\boldsymbol{p},\boldsymbol{Y}_{n+1}=\boldsymbol{Y}$, zero vectors $\boldsymbol{p}_{\text{new}}$ and $\mathbf{y}_{\text{new}}$ of size $T$.

10:     for $t=1$ to $T$ do

11:         Given $D_{n+1,t}=\{\boldsymbol{p}_{n+1},\boldsymbol{Y}_{n+1}\}$, draw a sample $p_{t}(\mathbf{x}_{n+1})$ from $Logitnormal(m(D_{n+1,t}),v(D_{n+1,t}))$, and then draw a sample $y_{t}(\mathbf{x}_{n+1})$ from a Bernoulli distribution with parameter $p_{t}(\mathbf{x}_{n+1})$.

12:         Update $\boldsymbol{p}_{n+1}=(\boldsymbol{p}_{n+1}^{\prime},p_{t}(\mathbf{x}_{n+1}))^{\prime}$, $\boldsymbol{Y}_{n+1}=(\boldsymbol{Y}_{n+1}^{\prime},y_{t}(\mathbf{x}_{n+1}))^{\prime}$, $(\boldsymbol{p}_{\text{new}})_{t}=p_{t}(\mathbf{x}_{n+1})$, and $(\mathbf{y}_{\text{new}})_{t}=y_{t}(\mathbf{x}_{n+1})$.      

13:     Set $\boldsymbol{p}_{\text{new}}^{(j)}=\boldsymbol{p}_{\text{new}}$ and $\mathbf{y}_{\text{new}}^{(j)}=\mathbf{y}_{\text{new}}$.

14:Take pointwise median from $\{\boldsymbol{p}_{\text{new}}^{(j)}\}_{j=1,\ldots,J}$ and $\{\mathbf{y}_{\text{new}}^{(j)}\}_{j=1,\ldots,J}$.