{'fraction_non_alphanumeric': 0.04103273397879207, 'fraction_numerical': 0.011331925942102836, 'mean_word_length': 4.345959268387599, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Latent function-on-scalar regression models for observed sequences of binary data: a restricted likelihood approach

In this paper, we study a functional regression setting where the random response curve is unobserved, and only its dichotomized version observed at a sequence of correlated binary data is available. We propose a practical computational framework for maximum likelihood analysis via the parameter expansion technique. Compared to existing methods, our proposal relies on the use of a complete data likelihood, with the advantage of being able to handle non-equally spaced and missing observations effectively. The proposed method is used in the Function-on-Scalar regression setting, with the latent response variable being a Gaussian random element taking values in a separable Hilbert space. Smooth estimations of functional regression coefficients and principal components are provided by introducing an adaptive MCEM algorithm that circumvents selecting the smoothing parameters. Finally, the performance of our novel method is demonstrated by various simulation studies and on a real case study. The proposed method is implemented in the R package dfrr.


Introduction
It is common that in some functional regression problems the response curve is observed as a sequence of correlated binary or multilevel data. This kind of situations can be handled via the family of generalized functional regression models. Bayesian methods are popular approaches for analyzing such data, see Goldsmith et al. [3], Meyer et al. [12] and van der Linde [18] among others. However, Bayesian methods show some limitations, including heavy computations and lengthy fitting procedures. They also often lack flexibility in implementing necessary constraints, such as the orthogonality of the eigen-functions, which is vital for the identifiability of the principal components. Typically, the likelihood analysis of such models is supposed not to be practicable, because of computational burden and difficulty of handling irregular and missing data [4]. Nonetheless, this paper aims at providing a likelihood analysis framework for probit functional regression models that is flexible enough to handle non-equally, irregular, and hence missing data. This work is motivated by the Madras-dataset [2, pp. 234-243], including a sequence of binary variables whose values indicate the presence or absence of Schizophrenia related symptoms in a set of patients observed during a time period. This dataset has been previously studied and analyzed by longitudinal methods [see 16]. Lately, the analysis of longitudinal data using methods from functional data analysis has been explored [7,21,23]. To get an insight into the relationship between functional and longitudinal data analysis we refer to Zhao et al. [22].
Alternative approaches to a likelihood analysis considered in literature include non-parametric methods [4], the use of approximating likelihood functions [19,17,8], or the use of simplified and relaxed assumptions [15]. Specifically, Hall et al. [4] proposed a non-parametric estimation of the mean and kernel function of the underlying Gaussian process for a correlated sequence of binary data. They argued that the maximum likelihood approach is computationally demanding and unstable, because it would require a large number of parameters to ensure a sufficiently flexible parameterization of the underlying Gaussian process. They also claim that the maximum likelihood approach is hard to use in the irregular data case, where its implementation requires imposing additional assumptions that limit its flexibility. However, the nonparametric method proposed in Hall et al. [4] depends on the very limiting assumption that the variation of the underlying process about its mean is relatively small, and assumption that is not needed in the maximum likelihood approach here proposed. Moreover, a major drawback of their work is that it does not handle covariates in modeling the location parameter.
Wang & Shi [19] proposed a generalized Gaussian process regression model for the binary case, in which location and the covariance operator are allowed depending on covariates. However, the approach is based on approximating the likelihood function by using a pre-specified kernel function that limits the flexibility of the model. A criticism of other likelihood-based approaches, like Scheipl et al. [15], is that their method ignores the possible correlation structure of the residual functions over the domain, which is typically present with functional data. Additionally, a reliable approach to the estimation of the principal components is not provided.
The estimation of the covariance operator, or equivalently of the functional principal components, is crucial in the analysis of longitudinal and functional data [6,9,13]. The functional principal components explain the pattern of variation of the observed random functions, and have many applications in dimension reduction [14]. Although often the covariance function is specified a priori [19,17], our likelihood analysis method does not rely on prior specification of the covariance function, which can be considered as a free parameter.
In conclusion, the method proposed in this paper has several advantages over the non-parametric, Bayesian, and other likelihood-based methods currently developed in the literature: it does not depend on additional unnecessary assumptions, it is computationally faster and more stable than Bayesian approaches, it does not ignore the correlation structure of the residual functions over the domain, and it also treats the covariance function as a free parameter.
In this paper, we consider a Function-on-Scalar Regression (FoSR) model setting, where the response function is a latent Gaussian process, and we only observe its realization as a sequence of correlated binary observations. We propose a maximum likelihood approach, inspired by the idea of parameter expansion [20], to carry out the estimation of regression coefficients and of the kernel function. Since we aim at automatically obtaining smooth functional parameters, functional regression coefficients and eigen-functions are estimated by introducing a novel method, namely the Adaptive Monte Carlo Expectation-Maximization (AMCEM) algorithm, which does not require the selection of the smoothing parameter. The AMCEM algorithm is a Monte Carlo EM algorithm in which the support of the distribution of the latent functions is restricted to a region consisting of functions that are smooth enough. It is adaptive in that the algorithm narrows the acceptance region step-by-step, so that the latent functions that fit the binary sequences become gradually smoother until the acceptance region converges to a non-empty set.
The rest of this paper is organized as follows. In Section 2, the dichotomized FoSR model is introduced and the identifiability issue of the parameters is thoroughly investigated. In Section 3, the estimation method and the AMCEM algorithm are described. The performance of the method is demonstrated by various simulation studies, and by showing the results of the analysis of the Madras dataset, in Section 4. Section 5 is devoted to discussion and conclusion. Proofs are collected in the Appendix.

The Functional Probit Regression Model
This article is concerned with functional data of the form Þ i ∈ L 2 [0, 1]; i = 1, 2, . . . , N . Let Þ i be a Gaussian variable with a mean function β x i and covariance operator T , where x i = (1, x i1 , . . . , x i(q−1) ) is a q-dimensional vector of covariates, and β = (β 0 , β 1 , . . . , β (q−1) ) is a q-dimensional vector of functions in L 2 [0, 1], representing the unknown regression coefficients. Suppose that the functions Þ i are latent and unobserved, that they are affected by measurement error, and only observed in a dichotomized and sparse version. That is, for each Þ i and some points 0 ≤ t i1 , t i2 , . . . , t iM i ≤ 1, observations are of the form: where I(·) is the indicator function. For describing the latent functions we thus consider a Function-on-Scalar Regression (FoSR) model of the form where The residual function ε i ∈ L 2 [0, 1] is a Gaussian variable with mean function zero and covariance operator T = j≥1 ν j ψ j ⊗ ψ j , where ⊗ denotes the tensor product and {ψ j ∈ L 2 [0, 1], j ∈ N} is a complete orthonormal basis of Image(T ). This means that T has kernel function K(s, t) = j≥1 ν j ψ j (s)ψ j (t). The error term it is independent of the residual function, and for each i = 1, 2, . . . , N and t ∈ [0, 1], it is independently distributed as a normal random variable with mean zero. The term it denotes the measurement error, while the residual function ε i (t) explains the complex and smooth structure of dependence along the domain within each sample. The variance of it together with the parameters β and T are not identifiable. The identifiability issues associated to the model parameters is discussed in Section 2.1.
The probability density function (pdf) of the binary response vector Y i = (Y i1 , Y i2 , . . . , Y iM i ) given the parameters β and T is . . .

Identifiability Issues
The parameters of the FoSR model (3) are not identifiable unless we add a constraint like K(t, t) = 1 for all t ∈ [0, 1]. In fact, consider the transformation is a linear operator and, for all g ∈ in which f is a positive and bounded function in , and covariance operator T = CT C having kernel function Since f is positive, by using the change of variable u = F i u in equation (4), we obtain Thus, (µ i , Σ i ), and hence (β, T, var( it )) are not identifiable. Although (β, T, var( it )) are not identifiable, a function of these parameters of the form (α, R) is identifiable, where α(t) = K(t, t) −1/2 β(t) and R = LT L. L is a linear mapping from L 2 [0, 1] to L 2 [0, 1] and, for any g ∈ L 2 [0, 1], we have (Lg)(t) = K(t, t) −1/2 g(t).
Remark 1. Note that α and R are respectively the mean function and covariance operator of the standardized version of the Gaussian process Þ in (3), given by and the kernel of the covariance operator R is given by Proposition 1. Consider the FoSR model given by (2)-(3), and let its functional response be dichotomized as in (1). Then, the parameters α and R are identifiable.
The relationship between the covariance operators of the standardized and non-standardized versions of the underlying Gaussian process is provided in the next proposition.
. Then, we have T = L RL and R = LT L with ||R|| ≤ 1.
For the sake of identifiability of the variance of the measurement error term, we consider the following standardized version of (2), We assume for simplicity that L it ∼ N (0, σ 2 ), independently for each i = 1, 2 . . . , N and t ∈ [0, 1]. This assumption is equivalent to that of the variance of it being proportional to the variance of Þ i at t in the unstandardized original model (2).
In Section 3.2, a truncated version of the model (7) as presented in (8) is used to provide Maximum Likelihood Estimation (MLE) of the location parameter, and Restricted Maximum Likelihood (REML) estimation of the covariance function. However, we first propose in Section 3.1 a framework related to the REML approach that facilitates providing unbiased estimation of covariance operators in infinite-dimensional Hilbert spaces.

Unbiased Estimation of the Covariance Operator
The REML approach is a special case of maximum likelihood estimation that uses a likelihood function provided from a transformed data, so that nuisance parameters have no effects on estimation. Although a likelihood function in the infinite-dimensional separable Hilbert spaces cannot be defined, an approach similar to REML can be used to provide an unbiased estimation of covariance function in such spaces.
First, note that some of the key results that hold for Gaussian random vectors in finite-dimensional spaces are also valid in the Gaussian case in infinite-dimensional separable Hilbert spaces.
Proposition 2. Let X 1 and X 2 be two jointly Gaussian functional variables in a separable Hilbert space H with mean functions m 1 and m 2 and covariance operators C 1 and C 2 , respectively. Then, X 1 and X 2 are independent if C 12 = C 21 = 0, where C 12 and C 21 are the cross-covariance operators.
Consider the matrix A = I N − X(X X) −1 X with rank(A) = N − q. Since A is not a full rank matrix, the eigen-value decomposition of matrix A is where U is a N × (N − q) matrix whose ith column is the eigen-vector u i of A, and Λ is the diagonal matrix Λ = diag(ν i ) 1≤i≤(N −q) whose diagonal elements are the corresponding positive eigen-values. Thus, noting that we have that Þ * = U Þ = U ε is a Gaussian random vector with mean vector of zero functions. From Proposition 3 (proof given in the appendix) below, it follows that the elements of Þ * are independent Gaussian functional variables with mean function zero and covariance operator T .

Proposition 3.
Suppose H is a separable Hilbert space and X = (X 1 , . . . , X n ) is a n-dimensional vector of H-valued independent Gaussian elements with mean functions zero and covariance operators C. let Q = (q 1 , q 2 , . . . , q m ) be a n × m matrix whose columns are orthogonal unit vectors, that is is a m-dimensional vector whose elements are independent Gaussian functional variables with mean functions zero and covariance operators C.
. From Proposition 3, Þ * k are N − q independent and identically distributed Gaussian random elements with mean function zero and covariance operators T . Hence, an unbiased estimate of T can be obtained by the N − q independent random samples Þ * where ⊗ stands for the tensor product 1 .
In the next section, we describe the EM algorithm, we provide the ML estimator of β, and we describe the REML estimator of T , similarly to what already proposed in equation (9).

Model inference via Expectation-Maximization
The likelihood function depends on the covariance function through the matrix Σ θ , thus the estimation of Θ and A and hence T are obtained by eigendecomposition of Σ θ . Consequently, to estimate the parameters of the FoSR model in (2), it is sufficient to estimate the Fourier coefficients of the functional parameters B, the covariance matrix Σ θ and σ 2 . Smooth estimations of β and ψ j , j ≥ 1 are given by introducing the AMCEM algorithm. In our estimation method, J in model (8) can be chosen large enough to reflect (7) with a high precision.
The EM algorithm is an iterative method to derive maximum likelihood estimates of parameters when the model depends on unobserved latent variables [1]. The maximization step of the EM algorithm concerning the standardized parameters does not yield a closed-form solution in the case of our model. Inspired by the parameter expansion technique proposed by Liu et al. [10], 1 Throughout the paper, we indicate with ⊗ both the tensor product and the Kronecker product, and we have so far explicitly mentioned which of the two was used. Hereafter, we will not explicitly specify the meaning of the notation ⊗, if it will be possible to infer it with no ambiguity from the mathematics. and Liu & Wu [11], we consider expanding the parameter space and estimating unstandardized parameters. The standardized parameters, which are uniquely identifiable, are then computed from the estimated unstandardized parameters.
Let e = (e 1 , e 2 , . . . , e J ) be the vector of the first J basis functions. Let (8), we can rewrite the models (3) and (2) in a truncated form as is the vector of independent measurement error terms, whose standardized version given by (6) is Gaussian distributed with mean vector zero and variance- Let us assume that the sequence of binary data for subject i, i = 1, . . . , N , . . , W N } be two sets of random vectors and Z = [Z 1 , Z 2 , . . . , Z N ] be a N × J random matrix, collecting the subject-specific i.i.d. model variables, and let y, w and z be the respective corresponding observations. Assuming that γ = (B, Σ θ , σ 2 ) is the model set of parameters, the complete-data log-likelihood function is given by where f (z, w | γ) denotes the joint pdf of Z and W given the parameters, and f (y | w, γ) is the pdf of Y given W and the parameter set.
The E-step at iteration m + 1 of the EM algorithm given the parameter estimate γ (m) at iteration m, involves evaluating the expectation where the latter expression is the expectation of the log-likelihood function given the observed variable Y. Considering that the idea behind our smoothing procedure for the functional parameters is to substitute Q(γ | y, γ (m) ) by ] 1≤j,k≤J is the penalty matrix associated with the basis functions e j (t), and e (n) j (t) is the n th derivative of e j (t). Throughout this paper, we set n = 2. The parameter λ controls the smoothness of the latent random functions, yielding to a smooth estimation of the eigen-functions.
Note that and thus Q λ (γ | y, γ (m) ) can be rewritten as where f (z | B, Σ θ ) stands for the pdf of Z and f (w | z, σ 2 ) is the pdf of W given Z.
To implement the M-step, we derive closed-form REML estimation of Σ θ and ML estimations of B and σ 2 . For the REML estimation of Σ θ , we substitute the restricted log-likelihood function to the log-likelihood function in (10). Differentiating with respect to Σ θ and setting the derivative equal to zero yields where Z * k = Z u k . Estimation of eigen-values and smoothed eigen-functions is then obtained by a simple eigen-decomposition of Σ ) are the eigen-values and eigen-vectors of Σ (m+1) θ , respectively.
By setting the first-order derivative of (10) with respect to B equal to zero, the update of the estimation of B is obtained as follows: Finally, the estimation of σ 2 at iteration m + 1 is updated by Let us define the matrix D as Then, the standardized regression coefficient α and the standardized covariance operator R can be updated by in which ϕ ) are respectively the eigenvalues and eigen-vectors of DΣ

The AMCEM Algorithm
In the E-step of the EM algorithm described above, it is necessary to compute the sufficient statistics in equations (11), (12) and (13), specifically: None of the above expectations have closed-form expressions, thus a Monte Carlo method needs to be employed to compute the integrations numerically. For the outer expectation, a Gibbs sampling algorithm can be developed to generate samples from the truncated multivariate normal distribution, and for the inner expectation, an accept-and-reject scheme can be used. The AM-CEM algorithm is "adaptive" because the sequence of smoothing parameters decreases at each step, inducing the latent functions fitted to the data to be smoother at each step proceeding along the algorithm. As a result, the AMCEM algorithm yields estimations of functional parameters which are as smooth as possible.
In order to apply the Gibbs sampling procedure, the full conditional distributions need to be specified. The conditional In equations (17) and (18) excluding the j th row and j th column. To generate samples from N (τ ij , σ 2 ij ) truncated to C ij , a simple method of inverse transform sampling can be employed. Let Φ x; τ ij , σ 2 ij be the cumulative distribution function of N (τ ij , σ 2 ij ), and define p = Φ 0; τ ij , σ 2 ij . Let u be a random sample from the uniform distribution U (0, p) if y ij = 0, and from U (p, 1) otherwise. A random sample from N (τ ij , σ 2 ij ) truncated to C ij is obtained by Φ −1 u; τ ij , σ 2 ij . To compute the inner expectations, consider that the conditional distribution of Z i given W i and γ (m) , for Thus, the inner expectation can be computed numerically by generating random samples from the given conditional distribution, truncated to the region specified by V λ . Random sample generation from this truncated multivariate normal distribution can be achieved by implementing an accept-and-reject scheme. The rejection rate of this procedure depends on λ i , with poor specification of λ i possibly leading to a high rejection rate. In the adaptive MCEM algorithm, λ i is chosen such that the rejection rate can be controlled at a prespecified level. The smoothing parameter λ i is decreasingly updated at each step of the Gibbs sampler. This approach avoids the selection of the smoothing parameter for the latent variables and eigen-functions, thus providing an automatically tuned smooth estimation.

Algorithm (AMCEM). Adaptive Monte Carlo EM
At iteration m + 1: and σ 2 (m+1) from the generated sample (w . Hence, the rejection rate of the algorithm can be controlled to be close to an a priori specified level δ. Note that, at the second iteration of the EM algorithm (see line 3), the parameters differ significantly from their initial values, and hence the distribution of Z i P (n) Z i | W i , γ (1) changes dramatically relative to the distribution of Z i P (n) Z i | W i , γ (0) at the first iteration, possibly causing the procedure at line 6 to yield sub-optimal results. We therefore choose to employ this procedure only from the third iteration of the algorithm.
− Line 14: obtaining smooth functions that fit well the sequence of binary data is guaranteed by the implemented accept-and-reject procedure.
The validation function in the AMCEM algorithm for a given 0 < δ < 1 is defined by andK i,δ is the evaluation of the covariance function estimated by the generated samples. Since when increasing the smoothing level the fitted latent functions cause the inflation of σ 2 , the validation function prevents the algorithm from over-smoothing.
The sequence λ (m) i is decreasing with respect to m, which causes the imputed missing functions to get smoother as the AMCEM proceeds. The decreasing sequence λ (m) i is bounded from below, thus implying its convergence. Furthermore, the convergence of λ (m) i corresponds to the convergence of the rejection rate δ (m) to zero. Since δ (m) converges to zero, V δ includes the support of Z given γ (m) and Y.
With the random samples generated at iteration m + 1, the parameters B, Σ θ and σ 2 are updated as follows: and u i is the i th column of U.

Further Comments
In the model description given in Section 3, the set of basis {e i } i≥1 was considered to be orthonormal. This assumption can be relaxed without any noticeable changes in our arguments, that is {e i } i≥1 can be chosen as any arbitrary basis, such as for example B-splines. The only differences are in estimating the standardized parameters in equations (14) and (15), and extracting the eigen-functions from the matrix Σ θ . The new standardized estimates of α and R are given by where Ω = [ e i , e j ] 1≤i,j≤J is a J × J matrix and ϕ ) are respectively the eigen-values and eigen-vectors of In our approach for estimating the functional parameters as described in Section (3.2), we proposed the smooth estimation of both functional regression coefficients and eigen-functions. As already mentioned, the smooth estimation of eigen-functions is obtained by introducing the novel AMCEM algorithm, which automatically tunes the smoothing parameters. It is worth noticing that the method, in addition to the eigen-functions, can also provide smooth estimations of the regression coefficients.
Prediction is often a crucial task in statistical modeling. Two different scenarios can be considered for prediction in our modeling framework. First, suppose that only the covariate values are given for a new individual. The estimation of the individual trajectory for the new case on its whole domain is then given byŴ =Ê[W | X = x] = (Be) x. Another possible prediction scenario entails that, in addition to the value of the covariates, the value of Y (t) is also observed for some time points t 1 , t 2 , . . . , t M . The individual trajectory for this new case can be then estimated bŷ The expectation E(W | X = x, Y = y) can be computed by MCMC integration approximation, using Gibbs sampling.

Applications
In this section, we evaluate the performance of the AMCEM algorithm via simulation studies and a real-world example. In the simulation studies, we examine the effect of: sample size, the number of points sampled per curve, the regularity/irregularity of the design, the magnitude of the variance of the measurement error, and the complexity of the covariance structure. For the regular/irregular designs, four scenarios are taken into account: (R) regular designs with equally spaced time points, (RT) regular designs with equally spaced right-truncated time points, (RM) regular designs with equally spaced and missing at random time points, (IRS) irregular designs with stochastic time points.

Simulation Studies
In our simulation study, we consider the FoSR model (2) with q = 1, given by in which ε ij are identically and independently distributed standard normal random variables. The parameters ν j and ψ j are eigen-values and eigenfunctions of the covariance operator of ε i , and also it iid ∼ N (0, σ 2 K (t, t)), where K(t, t) = j≥1 ν j ψ 2 j (t). A dichotomized version of this model is observed as y ij = I (W(t ij ) > 0) , i = 1, 2, . . . , N, j = 1, 2, . . . , M i .
In the data generating procedure, we consider two different sets of functional regression coefficients, detailed below. For what concerns the covariance structure of model (21), we assume that only the first p eigen-values are non-zero. The eigen-values are chosen as ν j = rρ j−1 for fixed r > 0 and 0 < ρ < 1, and the parameter ρ determines the decay rate of the eigen-values. If ρ tends to zero (one), the sequence of eigen-values decreases fast (slowly), thus leading to smoother (less smooth) residual functions, and to more (less) correlated binary variables in the sequence. For the eigen-functions, we set ψ 2j (t) = √ 2 sin (2πjt) and ψ 2j+1 (t) = √ 2 cos (2π (j + 1) t) for j ≥ 1. For what concerns the other model parameters, the structure of the simulation is as follows: • each of the four scenarios (R), (RT), (RM), and (IRS) is used as sampling design.
• N = 50, 100 is used respectively as small and large sample size.
• M = 12, 36 is used respectively as small and large number of sampling points per curve.
• σ 2 = 0.2, 0.8 is used respectively as small and large measurement error variance.
• ρ = 0.1, 0.4 is used respectively as large and small correlation within the binary variables in the sequence.
In the parameters estimation procedure, the minimal number of sampling points that is required to be measured per subject depends on the complexity of the functional parameters. Here we consider two cases for the functional regression coefficients: (1) a 'simple' case in which β 0 (t) = − cos(2πt) and β 1 (t) = 1 − 2t, and (2) a 'complex' case in which β 0 (t) = − cos(2πt) and β 1 (t) = − sin(4πt). The simulation study is run for all designs combinations, and for all combinations of the parameters N , M , σ 2 and ρ, by using an MCMC algorithm with 2,000 iterations. The accuracy of the AMCEM algorithm is evaluated by computing the mean square error (MSE) of estimations. The accuracy and efficiency are also compared to those of the alternative method implemented in the package pffr. Results are shown in Tables 1 and 2. To provide a better comparison of AMCEM and pffr, the unstandardized and standardized estimations of parameters are illustrated in Figures 1 and 4 for a single run of the simulation study. In the caption of each figure, a short combination of letters is given to indicate the combination of sampling design and complexity of regression coefficients the figure refers to. The first row of panels in each figure shows the unstandardized functional parameters, while the second row shows the standardized parameters. Since the unstandardized regression coefficients and the covariance operator are not identifiable, the functions illustrated in the first row are multiplied by a constant scalar value. The constant scalar is selected so that the unstandardized β 0 are similar in their L 2 -norm. For comparing the performance of AMCEM and pffr, we need to consider the standardized functional parameters, which are uniquely identifiable. To obtain a standardized version of pffr estimations, we need to first obtain an estimation of the covariance function, which is not provided directly by pffr. To fix this drawback, we used the residuals provided by pffr and the function fpca.sc in the package refund.
It can be noticed that in all simulation scenarios shown in the figures, e.g., scenario R(s), RT(s), R(c), and RT(c) with the prespecified aforementioned    simulation parameters, the unstandardized estimations provided by AMCEM are smoother than the ones obtained via pffr. Furthermore, the standardized estimations provided by AMCEM are more accurate. The results given in Table 1 and Table 2 show that the MSE of the estimations obtained both via AMCEM and via pffr decrease as the parameters N and M increase, as reasonably expected. Moreover, the precision of the estimation of the covariance operator (including the eigen-values and eigenfunctions) decreases when the value of σ 2 increases. This implies that the accuracy of AMCEM in estimating the covariance operator decreases (increases) as the total variation of the outcome is less (more) explained by the functional part of the model. By comparing Table 1 and Table 2, it can be noticed that an increase in ρ leads to a reduction in the MSE of the estimations of the functional Table 1: Mean square error of estimation of parameters (×10 4 ) with ρ = 0.1.  regression coefficients in the AMCEM method. As expected, if the value of the parameter ρ increases, the accuracy of pffr in estimating the functional regression coefficients increases dramatically. This can be explained by the fact that, as ρ approaches 1, the within-subject correlation of the sampling points decreases, an assumption on which the likelihood function in pffr model is based. In summary, according to the results presented in Table 1 and Table 2, AMCEM outperforms pffr in estimating the functional regression coefficients in a general situation. In connection to the results mentioned above, an increase in ρ causes a decrease in the accuracy of the estimation of the first principal component in AMCEM, and an increase in the accuracy when estimating the second. This can be explained by noting that an increase in ρ causes a decrease in the proportion of variation explained by the first eigen-function, and an increase in the proportion explained by the second. As a final note, we emphasize that the accuracy in the estimation of functional regression coefficients and functional principal components for incomplete designs is the highest in design (RT), followed by designs (RM) and (IRS), as expected from the pattern of missingness in the data.

Analysis of the Madras dataset
The Madras dataset (shortly, Madras-Data) is a collection of monthly records that indicates the presence or absence of a set of six positive (hallucinations, delusions, thought disorders) and negative (flat affect, apathy, withdrawal) psychiatric symptoms [2, pp. 234-243]. In this section, we apply our method to a dataset that consists of a sequence of binary observations indicating the presence or absence of 'thought disorder' in 86 patients over one year after initial hospitalization. The presence or absence of the symptom is recorded for each patient regularly every month, but for a few patients data is missing for some months. Each patient's record includes two features: age-at-onset and gender at the initial stage, which are considered as covariates in this study. The aim is to evaluate the effect of age and gender on the course of hospitalization. Concerning the response variable y ij , values of zero or one indicate presence or absence of the 'thought disorder' symptom respectively for the i th subject at month j, j = 0, 1, . . . , 11. The response y ij is missing for some random i and j, thus this dataset can be considered as data showing an irregular design, as described in Sec. 1. The collection of data is illustrated in Figure 5. It can be noticed that the sample size is not balanced in the four different groups defined by gender and age above/below 20 years old.
There are two scalar covariates in the dataset. The variable age is coded as zero and one, with 0 = age-at-onset ≥ 20 and 1 = age-at-onset < 20. The scalar variable gender is also coded as zero and one, with 0 = male, 1 = female. The major question regarding this study is to assess the effect of covariates on the course of illness. We apply a FoSR model to Madras-data, to investigate the relationship of age and gender to the 'thought disorder' symptom. In this model, the binary response Y indicates the presence or absence of the symptom, while the latent function Þ is a proxy for the thought disorder intensity (TDI). Our analysis also reveals the smooth principal components, which can provide valuable information about the latent pattern of variability of TDI.
To investigate the effect of Age and Gender on TDI, the following FoSR model is considered: This model is fitted to Madras-data by employing the AMCEM algorithm. Figure 6 shows the convergence of the estimation of the scalar-valued parameters of the model, including the measurement error σ 2 and the average degree of smoothness in the sample, defined as The left panel shows the estimated average degree of smoothness of the latent functions in the log scale along the AMCEM algorithm iterations, and the right panel shows the estimated measurement error. It can be noticed from this figure that the AMCEM algorithm converges. Furthermore, the measurement error is estimated asσ 2 = 0.2, which demonstrates that a major portion of the variation in the observed binary sequences is described by the functional part of the model. The strategy employed in the AMCEM algorithm for providing smooth estimations of functional parameters is based on fitting smooth functions to the sequence of binary data, consequently leading to smooth functional regression coefficients and covariance function estimates. The smooth estimations of the latent TDI functions fitted to four selected binary sequences is shown in Figure 7. The estimated smooth mean functions for the four different combinations of covariates are depicted in Figure 8. In this figure, (a) is the mean function of the females aged under 20, i.e, β 0 + β 1 + β 2 + β 3 , (b) is the mean function of the females aged 20 and over, i.e., β 0 + β 2 , (c) is the mean function of males aged under 20, i.e., β 0 + β 1 , and (d) is the mean function of males aged 20 and over, i.e., β 0 . By inspection of Figure 8, marked differences in the shapes of the four mean functions can be observed across all four groups. For instance, TDI is higher for younger females at the onset, although it is lower at the end of the course of illness; besides, the shape of the mean function is quite different in younger and older females, with the former showing a steady decline, and the latter instead a final increase in TDI at the end of the year. Moreover, the rate of reduction of TDI in older males is lower, meaning that TDI is less affected by hospitalization in this group. Another outstanding result that can be inferred from the estimated functional regression coefficients is that, while younger females experience a steady decrease in TDI, younger males show a quite different patter, with a stationary period followed by first an abrupt decrease in TDI in the middle of the year, and then by an increase in TDI towards the end.
The smooth standardized covariance function is illustrated in Figure 9, which carries valuable information on the illness trajectories. The first four smooth principal components of the covariance function are shown in Figure  10. The first two principal components account for 95.02% of the variation in TDI. These principal components reveal meaningful additional information for the illness trajectories. PC1 indicates that 73.07% of the variability of the trajectories is related to the variability of the TDI at baseline. PC2 indicates that 21.95% of the variability in the trajectories is related to the variability in the rate of TDI reduction. According to PC1 and PC2, one could consider a mixed-effect regression model with both random intercept and random slope instead of FoSR, to account for this pattern of variation. In order to compare the two different prediction scenarios described in Sec. 3.4, we consider four random samples with a few observed data. Predictions of the corresponding subject-specific trajectories given only the covariates, i.e., E(TDI | X = x) , and given the covariates and the observed sequence E(TDI | X = x, Y = y), are illustrated in Figure 11.

Conclusions and Discussion
We proposed a restricted maximum likelihood analysis framework for the FoSR model for the case in which a dichotomized version of the response curve is observed. A class of identifiable parameters was introduced, and a novel algorithm, namely the AMCEM algorithm, was proposed to provide smooth estimations of functional parameters. This approach does not rely on selecting hyperparameters, which is instead typical in penalized regression   Figure 11: prediction of the subject-specific latent functions for four not fully observed binary sequences. Dots in the plot indicate the sequence of observations of the 'thought disorder' symptom for the subject: empty circles for absence, filled circles for presence. Solid lines show the predicted latent functions, given the covariates only (grey) or given the covariates and the sequence of observed binary response (blue).
problems. In a quite extensive simulation study, we considered four different sampling designs including regular, truncated regular, regular with missing at random, and completely irregular designs for the observed timings of the observations. This simulation study demonstrated that the AMCEM algorithm provides satisfying results in all the aforementioned designs. We ran the simulation study for different combinations of sample sizes, number of points sampled per curve, magnitudes of the measurement error, and complexity of the covariance function. An R package named dfrr is available on CRAN for implementing our proposed method. In comparison to the pffr, the R package dfrr provides more accurate estimates; moreover, it also allows for the estimation of measurement error and covariance function, all giving useful information about the underlying process. Our proposed method was applied to the Madras longitudinal schizophrenia data to examine the effects of age and gender on the presence or absence of the 'thought disorder' symptom during the first year of hospitalization. Our analysis concludes that the major variation of the binary sequences can be described by a FoSR model. Comparing the estimated mean functions revealed that the younger females (older males) are more (less) affected by hospitalization with respect to the other groups. Another outstanding result was that the younger females experience a steady decrease in the TDI, while in younger males a quite different pattern is observed. By inspection of the estimated principal components one can understand the main pattern of functional variation in the data: 73.07% of the variation of the 'thought disorder' intensity is due to differences in the baseline, while 21.95% of the variation is due to differences in the rate of decrease of the 'thought disorder' intensity during hospitalization.
The current work can be extended in two directions. The covariates on the right-hand side of the model can be extended to functional covariates, and the left-hand side of the model can either be any discretized version of the functional response or censored from below or above. Thus, ordinal and Tobit function-on-function regression models are the natural extensions of this work. and for any g ∈ L 2 [0, 1] (L K CT CL K g) (t) = K (t, t) − 1 2 f (t) K(s, t)f (s)K (s, s) − 1 2 g(s)ds = f (t) −1 f (t)K(t, t) − 1 2 K(s, t)f (s)f (s) −1 K(s, s) − 1 2 g(s)ds = K(t, t) − 1 2 K(s, t)K(s, s) − 1 2 g(s)ds = (Rg) (t).
Proof. Proposition 2. Let X 1 and X 2 be independent variables. Suppose f and g are two functions in H and µ X 1 ,X 2 is a borel probability measure on the Hilbert space H 2 = H × H, then Therefore, C 12 f, g = 0 which implies C 12 f = 0 for all f ∈ H. With regard to C 12 f = 0, we have C 12 = 0. Similarly, we can prove C 21 is equal to zero. Conversly, we assume C 12 = 0. Thus, for all f, g ∈ H, we have C 12 f, g = 0 which shows cov( X 1 − m 1 , f , X 2 − m 2 , g ) = E X 1 − m 1 , f X 2 − m 2 , g = 0. (26) Noting that Y 1 = X 1 − m 1 , f and Y 2 = X 2 − m 2 , g are jointly Gaussian, from (26), Y 1 and Y 2 are independent. We can write As a result, X 1 and X 2 are independent.
We know that Y i is Gaussian with the mean function zero. Thus, if f, g ∈ H, we can write n l=1 q ik q jl E X k , f X l , g = n k=1 q ik q jk E X k , f X k , g = Cf, g q i q j = Cf, g δ ij , which shows C ii = C and C ij are equal to zero, for all i = j, where C ij is the cross-covariance operator of Y i and Y j . As a result, based on Proposition 2, Y i for i = 1, 2, . . . , m are independent.