To derive the GFPCA solution,
we first center the observed \(Y_i(t)\)  based on their marginal
mean \(\mu _Y(t) = \mathsf {E}[Y_i(t)]\) ,
estimated through a simple smoother of all \((t_{ij}, Y_i(t_{ij}))\) -pairs via a
generalized additive model [1]} with
response function \(g(\cdot )\)  and the appropriate exponential family for the response.
The covariance of the latent process can then be approximated by
\(\widehat{\left[ X_i(s), X_i(t) \right] \approx \frac{\hat{\sigma }_Y(s,t)}{g^{(1)}[\mu _X(s)] \cdot g^{(1)}[\mu _X(t)]},}with \sigma _Y(s,t) = \mathsf {E}[Y_{c,i}(s) \cdot Y_{c,i}(t)] based on the centered curves Y_{c,i}(t),the marginal mean \mu _X(t) estimatedaccordingly to \mu _Y(t) and g^{(1)}(\cdot ) the first derivative of the response function.For given time points s_1 and s_2, \sigma _Y(s_1,s_2)is estimated as the mean of all pairwise products y_{c,i}(s_1)\cdot y_{c,i}(s_2).The estimated surface \hat{\sigma }_Y(s,t) is a smoothed version of \sigma _Y(s,t) using a bivariate tensor product P-spline basis \cite {fahrmeir_2013}. Since the surface is expected to show some discontinuity along the diagonalthis smoothing step is performed under exclusion of the diagonal elements\cite {yao_2005}.The FPCs \psi _k(t) and their respective eigenvalues \tau _kare then estimated from the spectral decomposition of \widehat{\left[X_i(t), X_i(s)\right].Our approach deviates slightly from the method of \cite {hall_2008} it is based on: We mean-center the data beforetaking their crossproducts instead of subtracting the crossproduct of the estimated mean from the crossproducts of the data.In our experience, this yields smoother estimates of the covariance surface which are more amenable to a low-rank FPC representation.}\) 
