Thinking. 1. **Identify the Goal:** * Target Algorithm: `sae::pbmseSTFH` from the `sae` R package. * Target Language: Rcpp (C++ integrated with R). * Purpose: Estimate mean squared error (MSE) of the spatio-temporal Fay-Herriot (STFH) model using parametric bootstrap (PB). * Requirement: More efficient and faster than the existing R implementation. 2. **Analyze the Target: `sae::pbmseSTFH`** * What does it do? Fits a spatio-temporal Fay-Herriot model, then runs a parametric bootstrap to compute MSE estimates. * Model: STFH is a mixed model. Usually state-space or area-level. * Let's recall the exact model and estimation procedure from the `sae` package (or the underlying theory, Battese, Harter, Fuller? No, STFH is standard spatio-temporal small area estimation). * First step: Estimate model parameters (variance components) using REML. * Second step: Calculate EBLUP (Empirical Best Linear Unbiased Predictor). * Third step: Parametric bootstrap to estimate MSE. * Let's check the `sae::pbmseSTFH` function signature to understand inputs. * Formula for fixed effects. * Domain index, time index. * Variance of direct estimates. * Number of bootstrap replicates B. * Let's look at the structure of `sae::pbmseSTFH` (from memory or logical deduction, as I can't execute code). * Inputs: `formula`, `vardir`, `D` (number of domains), `T` (number of time periods), `B` (bootstrap replicates). * Output: `mse` (MSE estimates), `est` (EBLUP). * Core computational burden in R implementation: 1. REML optimization (maximizing profile likelihood) - usually done via `nlminb` or `optim`. 2. Computing EBLUPs (matrix algebra, inverses). 3. Bootstrap loop `for(b in 1:B)`: Generate new data from fitted model, re-estimate parameters (REML), compute new EBLUPs, compute bootstrap MSE. * Speed bottlenecks in the R code: * The REML optimization in the bootstrap loop. This is the heaviest part. * Large matrix inversions and multiplications within EBLUP computation. * R loops are slow compared to C++ loops. * Wait, `sae::pbmseSTFH` is already written in R. Let's verify the function if possible. (It's primarily R code, using some C, maybe, but the main loops and optimization are R). 3. **Deep Dive into the Algorithm (Standard STFH model)** * Model: $Y_{it} = X_{it}^T \beta + u_i + v_{it} + \epsilon_{it}$ * $u_i \sim N(0, \sigma_1^2)$ (domain-specific random effect) * $v_{it} \sim N(0, \sigma_2^2)$ (time-specific random effect, or domain-time interaction, or AR(1) etc. Let's check `pbmseSTFH` specification). * In `sae`, the STFH model is: $Y_{it} = X_{it}^T \beta + u_i + v_{it} + e_{it}$, where $v_{it} = \rho v_{i,t-1} + \eta_{it}$ (AR1). This is the standard "spatio-temporal" FH model. Let's stick to the model *exactly* as implemented in `sae::pbmseSTFH`. * Actually, `sae` package STFH model (from Molina, I., Saei, A. - maybe not, actually it's from Marhuenda et al.?) Let's trace the typical STFH model. * "STFH" in `sae` (vignette, help file): * `pbmseSTFH(formula, D, T, vardir, B = 200, ...)` * Model: $\theta_{it} = x_{it}^T \beta + u_i + v_{it} + \epsilon_{it}$. * $u_i \sim N(0, \sigma_{u}^2)$ (domain effect). * $v_{it}$ is either AR1 ($v_{it} = \rho v_{i,t-1} + \eta_{it}$, $\eta_{it} \sim N(0, \sigma_{v}^2)$) or random walk, or independent. Let's look at the function description. * Ah, `STFH` is "Spatio-Temporal Fay-Herriot" (Marhuenda, Molina, Morales, 2013). * Model formulation in `sae`: * Sampling model: $\hat{\theta}_{it} = \theta_{it} + e_{it}$, $e_{it} \sim N(0, \psi_{it})$. * Linking model: $\theta_{it} = x_{it}^T \beta + v_{it}$, where $v_{it} = \rho v_{i,t-1} + u_{it}$, with $u_i \sim N(0, \sigma_{u}^2)$ and $v_{i1} \sim N(0, \sigma_{v}^2/(1-\rho^2))$. * Wait, actually the standard ST model is: * $\theta_{it} = x_{it}^T \beta + u_i + v_{it}$ (domain + domain-time) * Where $v_{it} = \rho v_{i,t-1} + \epsilon_{it}$ (AR1). * The variance is $\Sigma_i = \sigma_u^2 J_T + \sigma_v^2 \Omega_i(\rho) + diag(\psi_i)$. * Often, $u_i$ is domain effect, $v_{it}$ is an AR1 process. * Optimization: REML / ML. Nelder-Mead or L-BFGS-B? * Let's check the actual `stEblup` or `pbmseSTFH` from the `sae` package. * Let's reconstruct the algorithm structure. 1. Build matrices. 2. Estimate $\hat{\beta}$, $\hat{\sigma}_u^2$, $\hat{\rho}$, $\hat{\sigma}_v^2$ (REML). 3. Calculate EBLUP. 4. For b in 1:B: * Generate $u_i^{(b)} \sim N(0, \hat{\sigma}_u^2)$, $v_{it}^{(b)}$ from the AR1 process with $\hat{\rho}$, $\hat{\sigma}_v^2$. * Generate $\theta_{it}^{(b)} = x_{it}^T \hat{\beta} + u_i^{(b)} + v_{it}^{(b)}$. * Generate $\hat{\theta}_{it}^{(b)} = \theta_{it}^{(b)} + e_{it}^{(b)}$, where $e_{it}^{(b)} \sim N(0, \psi_{it})$. * Calculate bootstrap estimators $\hat{\beta}^{(b)}$, $\hat{\sigma}_u^{2(b)}$, $\hat{\rho}^{(b)}$, $\hat{\sigma}_v^{2(b)}$ using REML (or ML) from the bootstrap data. * Calculate $\hat{\theta}_{it}^{(b)}$ (bootstrap EBLUP). * Compute bootstrap MSE: $(\hat{\theta}_{it}^{(b)} - \theta_{it}^{(b)})^2$. 5. Average over B. 4. **Feasibility and Efficiency Gains in Rcpp** * The bootstrap loop and the REML optimization are the main computational targets. * Steps in R: 1. **REML optimization**: Maximizing the (restricted) log-likelihood. * Function evaluation involves computing $\Sigma_i = \sigma_u^2 J_T + \sigma_v^2 \Omega_i(\rho) + diag(\psi_i)$. * Computing the determinant $\log|\Sigma_i|$. * Computing the quadratic form $(y_i - X_i\beta)^T \Sigma_i^{-1} (y_i - X_i\beta)$. * Calculating Profile ML/REML. * Fisher scoring / Newton-Raphson / BFGS / Nelder-Mead. * In R, `optim` or `nlminb` is called. * We can implement the REML optimization in C++ using the `nlopt` library (via `RcppNLopt`), `RcppEigen`/`RcppArmadillo` for the matrix algebra, or even implement a custom Fisher scoring algorithm. * Wait, could we just use `RcppEigen`/`RcppArmadillo` to compute the profile likelihood function and then let R's optimizer call it? * Yes, this is usually the best approach. Wrap the profile log-likelihood in a `fun` function in C++, compile it to an R function via `Rcpp::cppFunction` or `sourceCpp`. * Then R calls `optim(fn = cpp_fun, ...)`. * Can we do better? Remove the R overhead from the optimization loop entirely. * Can we write the whole `pbmseSTFH` in Rcpp? * Yes. * Algorithm: 1. Input: `y`, `X`, `D`, `T`, `vardir`, `B`. 2. `Eigen::MatrixXd` or `arma::mat` for data. 3. Function: `STFH_REML(theta, y, X, vardir, D, T)`. * Let $\sigma_u^2 = \exp(\theta_0)$, $\sigma_v^2 = \exp(\theta_1)$, $\rho = \text{tanh}(\theta_2)$ or transform appropriately. * Wait, how does R's `sae::pbmseSTFH` parameterize? It uses $\sigma_u^2$, $\sigma_v^2$, $\rho$. * Compute $\Sigma_i$ for each domain $i$. * Since $\Sigma_i$ is $T \times T$, $D$ matrices. * Cholesky decomposition of $\Sigma_i$: $L_i L_i^T = \Sigma_i$. * Determinant: $2 \sum \log L_{ii}$. * Inverse: $\Sigma_i^{-1} y_i = L_i^{-T} L_i^{-1} y_i$. * Calculate $\hat{\beta}_{GLS} = (X^T \Sigma^{-1} X)^{-1} X^T \Sigma^{-1} y$. * Calculate REML. 4. Optimize using `nlopt` (`LN_BOBYQA` or `LN_NELDERMEAD`). 5. Calculate EBLUP. * $\hat{\theta}_i = X_i \hat{\beta} + G_i Z_i^T \Sigma_i^{-1} (y_i - X_i \hat{\beta})$. * Where $G_i$ is the covariance matrix of the random effects (domain effect + AR1). This is a $T \times T$ matrix. * $Z_i$ is design matrix for random effects (identity for $D$ domains + T time? No, STFH has domain effect and domain-time AR1). * EBLUP calculation: $\hat{\theta} = X \hat{\beta} + \Sigma_v Z^T V^{-1} (y - X\hat{\beta})$. 6. Bootstrap: * For b = 1 to B: * Generate $u_i \sim N(0, \hat{\sigma}_u^2)$. * Generate $v_{i1} \sim N(0, \hat{\sigma}_v^2 / (1-\hat{\rho}^2))$. * Generate $v_{it} = \hat{\rho} v_{i,t-1} + \eta_{it}$, $\eta_{it} \sim N(0, \hat{\sigma}_v^2)$. * $\theta^b = X \hat{\beta} + u + v$. * $y^b = \theta^b + e$, $e \sim N(0, vardir)$. * Run the REML estimation on bootstrap data $(y^b, X, vardir)$. * Calculate bootstrap EBLUP $\hat{\theta}^b$. * Accumulate $(\hat{\theta}^b - \theta^b)^2$. 7. Return MSE estimates. 5. **Crucial Optimizations in C++** * **Parallelization**: The bootstrap loop is an obvious target for parallelization (OpenMP via `RcppParallel` or `Rcpp` with OpenMP). Since each bootstrap replicate is independent, we can parallelize the REML estimation and EBLUP computation for each replicate. * **REML Speed**: The bootstrap REML is the bottleneck. We can use starting values from the original fit to speed up convergence. * **Avoiding R memory management**: R's garbage collector can slow things down if we constantly call back. Pure C++ loop with Eigen/Armadillo is much faster. * **Profile likelihood function re-uses**: * $X$ is constant. * $vardir$ is constant. * Structure of $\Sigma_i$ is known. * We can exploit the block diagonal structure of $\Sigma$ (especially if $\sigma_u^2$ and $\sigma_v^2$ are globally estimated). * In standard FH, $\Sigma = diag(\Sigma_1, ..., \Sigma_D)$. * The REML function: * $l_{REML}(\sigma_u^2, \sigma_v^2, \rho) = const - \frac{1}{2} \sum_{i=1}^D \log |\Sigma_i| - \frac{1}{2} \log |X^T \Sigma^{-1} X| - \frac{1}{2} y^T P y$. * Where $P = \Sigma^{-1} - \Sigma^{-1} X (X^T \Sigma^{-1} X)^{-1} X^T \Sigma^{-1}$. * Computing this efficiently is key. * We can create a class that pre-allocates matrices and uses efficient Eigen routines. * **Matrix Inversion Speed**: * $\Sigma_i$ is $T \times T$. $T$ is usually small (like 10-20). So inversion is $O(T^3)$. * Number of inversions per REML evaluation: $D$. * Total per bootstrap: $D * N_{iter}$. * This is fine in C++ with Eigen/Armadillo compared to R which scales like $O(D T^3 * N_{iter})$ but with R overhead. * Vectorization: We can stack the data, but the covariance structure is different per domain if $\psi_{it}$ varies. The correlation structure from $u_i$ and $v_{it}$ is the same across domains. $\Sigma_i = \sigma_u^2 J_T + \sigma_v^2 \Omega + diag(\psi_i)$. So $D$ different $\Sigma_i$. 6. **Detailed Algorithm Translation** * **Inputs**: `y` (vector, length $n = D \times T$), `X` (matrix, $n \times p$), `D` (int), `T` (int), `vardir` (vector, length $n$), `B` (int). * **Helper Functions**: * `make_Sigma(sig_u2, sig_v2, rho, vardir_i)`: Returns $\Sigma_i$ (T x T). * `log_det_Sigma`: Sum of log diag of Cholesky. * `log_det_XtSinvX`: Using QR/chol of $X^T \Sigma^{-1} X$. * `neg2_log_REML(params, y, X, vardir, D, T)`: Return value for optimzer. * `ONE_REML_ESTIMATION(y, X, vardir, D, T, starting_vals)`: Perform optimization (e.g. BOBYQA from NLopt) and return $(\hat{\beta}, \hat{\sigma}_u^2, \hat{\sigma}_v^2, \hat{\rho}, \hat{\theta})$. * **Main Function**: * Estimate parameters from original data: `params0 = ONE_REML_ESTIMATION(y_orig, X, vardir, D, T)`. * Calculate $\hat{\theta}_{orig}$ EBLUP. * Parallel `for` loop over `b` (using OpenMP, e.g. `#pragma omp parallel for`): * Generate $\theta^b$, $y^b$. * `params_b = ONE_REML_ESTIMATION(y_b, X, vardir, D, T, params0)` (start from original params for speed). * Calculate $\hat{\theta}^b$. * Accumulate squared errors $(\hat{\theta}^b_i - \theta^b_i)^2$ in a shared vector `mse_sum`, protected by `#pragma omp atomic update`. * Return `mse = mse_sum / B`. 7. **Refining the REML Optimization in C++** * Using `nlopt` (`nloptr` R package or `RcppNLopt`). * `nlopt` solves minimization problems. * We minimize $-2 \log L_{REML}$. * Parameters need bounds: $\sigma_u^2 >= 0$, $\sigma_v^2 >= 0$, $\rho \in (-1, 1)$. * Using log-parameterization to avoid boundary issues and improve convergence: * $\theta_1 = \log(\sigma_u^2)$ * $\theta_2 = \log(\sigma_v^2)$ * $\theta_3 = \text{atanh}(\rho)$ * Wait, what does the original `sae::pbmseSTFH` use? REML, `nlminb` or `optim`. * Implement a Newton or Fisher scoring? The analytical gradient of REML for STFH is available (see Rao & Molina). This would allow much faster optimization! * Gradient of REML: * $\frac{\partial l_{REML}}{\partial \sigma_u^2} = -\frac{1}{2} \sum_i tr(P_i \frac{\partial \Sigma_i}{\partial \sigma_u^2})$ (or similar). * Let's derive or lookup the marginal REML for STFH. * Actually, the characteristic of the STFH model is that it can be written as $\hat{\theta}_{it} = x_{it}^T \beta + u_i + v_{it} + e_{it}$. * $V_i = \sigma_u^2 1_T 1_T^T + \sigma_v^2 \Omega_i + D_i$. * REML: $l_{REML} = const - \frac{1}{2} \sum_i \log|V_i| - \frac{1}{2} \log |\sum_i X_i^T V_i^{-1} X_i| - \frac{1}{2} \sum_i (y_i - X_i \hat{\beta})^T V_i^{-1} (y_i - X_i \hat{\beta})$. * Fisher scoring (FS) is often used in practice in SAE packages (like `sae`, `emdi`, `msae`). FS requires the score vector and expected Fisher information matrix. * Score vector elements: * $S(\sigma_u^2) = -\frac{1}{2} \sum_i tr(P_i \frac{\partial V_i}{\partial \sigma_u^2}) + \frac{1}{2} (y - X\hat{\beta})^T P_i \frac{\partial V_i}{\partial \sigma_u^2} P_i (y - X\hat{\beta})$ (No, the second term is zero if we use REML, wait. The score for REML is $S(\phi) = -\frac{1}{2} tr(P \frac{\partial V}{\partial \phi}) + \frac{1}{2} y^T P \frac{\partial V}{\partial \phi} P y$. * $\frac{\partial V_i}{\partial \sigma_u^2} = J_T$. * $\frac{\partial V_i}{\partial \sigma_v^2} = \Omega_i$. * $\frac{\partial V_i}{\partial \rho} = \sigma_v^2 \frac{\partial \Omega_i}{\partial \rho}$. * This is very standard. We can implement FS in C++ with Eigen. It will be extremely fast. * Let's just use a derivative-free optimizer (BOBYQA) for simplicity, OR FS. * If we use FS: * Initialize $\phi^{(0)}$. * For k=1, ...: * Compute $V_i$, $V_i^{-1}$. * Compute $\hat{\beta} = (\sum_i X_i^T V_i^{-1} X_i)^{-1} \sum_i X_i^T V_i^{-1} y_i$. * Compute $P_i = V_i^{-1} - V_i^{-1} X_i (\sum_i X_i^T V_i^{-1} X_i)^{-1} X_i^T V_i^{-1}$. * Compute Score $S(\phi)$. * Compute Fisher Info $\mathcal{I}(\phi)$. * $\phi^{(k+1)} = \phi^{(k)} + \mathcal{I}^{-1} S(\phi^{(k)})$. * Eigen implementation of FS will be incredibly tight. * Let's compare FS vs NLOPT. * FS needs gradient and expected information. Code is medium complexity. * NLOPT just needs the function value. Code is simple. * FS might require fewer iterations and be very fast. Let's do FS! * Wait, ensuring the algorithm converges (especially with many replications) is crucial. FS can overshoot (step too large). We can dampen the FS step (score weighting, line search). * Let's do the NLopt version first, as it's more robust and guarantees convergence for tricky datasets without careful Hessian modifications. But I can write an FS version too. * Actually, the `joineRML` package and others use FS. It usually just works if we truncate the step. We can also just use L-BFGS from NLopt. L-BFGS is very good and requires only gradient! * Let's provide the gradient of the REML to L-BFGS. This gives robustness and speed. 8. **Code Structure** * `RcppSTFH` class or just a set of C++ functions within `sourceCpp`. * Let's write the answer as a standalone `Rcpp` function file. * Functions: 1. `sig_u2_from_theta(theta[0])`, `sig_v2_from_theta(theta[1])`, `rho_from_theta(theta[2])`. 2. `mat Sigma_i(sig_u2, sig_v2, rho, vardir_i, T)`. 3. `double neg2_log_REML(const Eigen::VectorXd& params, ...)`. 4. `Eigen::VectorXd grad_REML(const Eigen::VectorXd& params, ...)`. 5. `Rcpp::List one_bootstrap(...)`. 6. `// [[Rcpp::export]]` `Rcpp::List pbmseSTFH_fast(Rcpp::NumericVector y, Rcpp::NumericMatrix X, int D, int T, Rcpp::NumericVector vardir, int B)`. * Let's detail the gradient of REML for L-BFGS. * $\phi = (\sigma_u^2, \sigma_v^2, \rho)$ or better $\phi = (\log \sigma_u^2, \log \sigma_v^2, \text{atanh}(\rho))$. * $V_i = \sigma_u^2 J_T + \sigma_v^2 \Omega_i + \text{diag}(\psi_i)$. * Let $\phi = (\phi_1, \phi_2, \phi_3)$. * $\sigma_u^2 = \exp(\phi_1)$, $\frac{\partial V_i}{\partial \phi_1} = \frac{\partial V_i}{\partial \sigma_u^2} \cdot \sigma_u^2 = J_T \cdot \sigma_u^2$. * $\sigma_v^2 = \exp(\phi_2)$, $\frac{\partial V_i}{\partial \phi_2} = \frac{\partial V_i}{\partial \sigma_v^2} \cdot \sigma_v^2 = \Omega_i \cdot \sigma_v^2$. * $\rho = \frac{\exp(2\phi_3) - 1}{\exp(2\phi_3) + 1}$. Let's use $\rho = \tanh(\phi_3)$. $\frac{\partial \rho}{\partial \phi_3} = \text{sech}^2(\phi_3) = 1 - \tanh^2(\phi_3) = 1 - \rho^2$. * $\frac{\partial V_i}{\partial \phi_3} = \frac{\partial V_i}{\partial \rho} \cdot \frac{\partial \rho}{\partial \phi_3} = \sigma_v^2 \frac{\partial \Omega_i}{\partial \rho} \cdot (1 - \rho^2)$. * The AR(1) matrix $\Omega_i$: * $(\Omega_i)_{st} = \frac{\rho^{|s-t|}}{1 - \rho^2}$. * Wait, in standard STFH, the AR(1) process is $v_{it} = \rho v_{i,t-1} + \eta_{it}$, $\eta_{it} \sim N(0, \sigma_v^2)$. * If $\text{Var}(v_{i1}) = \sigma_v^2 / (1 - \rho^2)$, then $\text{Var}(v_{it}) = \sigma_v^2 / (1 - \rho^2)$ for all t. * $\text{Cov}(v_{it}, v_{is}) = \frac{\sigma_v^2}{1 - \rho^2} \rho^{|t-s|}$. * Let $\Gamma_i = \sigma_v^2 \Omega_i$, where $\Omega_i = \frac{1}{1 - \rho^2} \begin{pmatrix} 1 & \rho & \cdots \\ \rho & 1 & \cdots \\ \vdots & \vdots & \ddots \end{pmatrix}$. * Actually, the model $v_{it} = \alpha v_{i,t-1} + u_{it}$ implies $\text{Cov}(v) = \sigma_u^2 / (1 - \alpha^2) \Gamma$ where $\Gamma_{st} = \alpha^{|s-t|}$. * So $V_i = \sigma_u^2 J_T + \frac{\sigma_v^2}{1 - \rho^2} C(\rho) + \text{diag}(\psi_i)$, where $C(\rho)_{st} = \rho^{|s-t|}$. * Let's check the exact model in `sae`: * "domain-time random effects $v_{it}$ follows an AR(1) process". * In `pbmseSTFH`, the model is $\theta_{it} = x_{it}^T \beta + u_i + v_{it}$, where $v_{it} = \phi v_{i,t-1} + \epsilon_{it}^{(2)}$. * The sampling error is $e_{it}$. * Variances: $var(u_i) = \sigma_1^2$, $var(\epsilon_{it}^{(2)}) = \sigma_2^2$. * So $V_i = \sigma_1^2 J_T + \frac{\sigma_2^2}{1-\phi^2} \tilde{\Omega}(\phi)$, where $\tilde{\Omega}_{st} = \phi^{|s-t|}$. * Let $\rho = \phi$. * Yes, this matches the standard parameterization. * Let's derive $\frac{\partial \Omega_i}{\partial \rho}$. * $\Omega_i = \frac{1}{1-\rho^2} C(\rho)$ where $C_{st}(\rho) = \rho^{|s-t|}$. * $\frac{\partial \Omega_i}{\partial \rho} = \frac{2\rho}{(1-\rho^2)^2} C(\rho) + \frac{1}{1-\rho^2} \frac{\partial C(\rho)}{\partial \rho}$. * $\frac{\partial C(\rho)_{st}}{\partial \rho} = |s-t| \rho^{|s-t| - 1}$ (if $|s-t| > 0$), else 0. * This is computable. We can precompute $C$ and $C'$. * OK, let's simplify the code by using a derivative-free or gradient-based optimizer available in RcppNLopt. * Actually, writing a complete, robust, efficient implementation from scratch requires a lot of matrix algebra. Let's make sure the code is solid. * **Alternative: RcppEigen using Cholesky updates for the bootstrap?** * In the bootstrap, the sampling error $e_{it}$ is generated. The $\psi_{it}$ is known and fixed. * The REML is run on $y^b$. * We can speed up the bootstrap REML by using starting values from the original dataset. * Can we speed it up more? The $X$ matrix is constant! The $\psi_i$ vector is constant! * Wait, if $\psi_i$ is constant, and we use the log-REML, the inverse $V_i^{-1}$ depends on $y^b$ only through the quadratic forms. * We can compute $V_i^{-1}$ and $V_i^{-1} X_i$ once per REML iteration. * Let's just use Newton's method (Fisher Scoring) since it converges quadratically and we can parallelize the score vector / information matrix computation. * Derivation of Fisher Scoring for STFH (REML): * $l_{REML} = -\frac{1}{2} \log |V| - \frac{1}{2} \log |X^T V^{-1} X| - \frac{1}{2} y^T P y$. * $K = X^T V^{-1} X$. * $\hat{\beta} = K^{-1} X^T V^{-1} y$. * $P = V^{-1} - V^{-1} X K^{-1} X^T V^{-1}$. * Score for $\sigma_k^2$ (general variance component): * $U(\sigma_k^2) = -\frac{1}{2} tr(P V_k) + \frac{1}{2} y^T P V_k P y$, where $V_k = \frac{\partial V}{\partial \sigma_k^2}$. * Wait, $V_k = \frac{\partial V}{\partial \phi_k}$. * Expected Information (Fisher Information): * $\mathcal{I}(\sigma_k^2, \sigma_l^2) = \frac{1}{2} tr(P V_k P V_l)$. * Since $V$ is block diagonal over domains: * $V = \text{diag}(V_1, \dots, V_D)$. * $U(\phi) = \sum_{i=1}^D U_i(\phi)$. * $\mathcal{I}(\phi)_{kl} = \sum_{i=1}^D \mathcal{I}_{i, kl}(\phi)$. * $U_i(\phi_k) = -\frac{1}{2} tr(P_i V_{i,k}) + \frac{1}{2} e_i^T P_i V_{i,k} P_i e_i$, where $e_i = y_i - X_i \hat{\beta}$. * $\mathcal{I}_{i, kl}(\phi) = \frac{1}{2} tr(P_i V_{i,k} P_i V_{i,l})$. * This is vectorizable perfectly. * Let's commit to Fisher Scoring in C++. * It avoids optimizer overhead from NLopt. * It is analytically correct and fast for well-behaved models. * Standard in `sae` package (check the source code, it uses Fisher scoring). * Wait, let's check if `sae::pbmseSTFH` uses Fisher Scoring. YES! The help page for `pbmseSTFH` says: * "REML is carried out by Fisher-scoring algorithm". * Ah! In that case, we just port it to C++! Excellent. * **Algorithm from `sae::pbmseSTFH.source` (reverse eng/remembering)**: * Function `stFH` or similar. * Parameters: $\sigma_u^2$, $\sigma_v^2$, $\phi$ (autocorrelation). * Iterates until convergence of variance components. * Let's write this cleanly in C++. * Steps in C++ Fisher Scoring Loop: 1. Calculate $\rho = \phi$. Initialize $\phi^0$. 2. For each iteration: * Compute $V_i$ for all $i$. * Cholesky: $L_{V_i}$ lower. * Compute $V_i^{-1} y_i$. * Compute $V_i^{-1} X_i$. * Accumulate $X^T V^{-1} y$, $X^T V^{-1} X$. * Compute $\hat{\beta} = (X^T V^{-1} X)^{-1} (X^T V^{-1} y)$. * Compute $e_i = y_i - X_i \hat{\beta}$. * Compute $V_i^{-1} e_i$ (we already have $V_i^{-1} y_i$ and $V_i^{-1} X_i$, so $V_i^{-1} e_i = V_i^{-1} y_i - V_i^{-1} X_i \hat{\beta}$). * Compute $P_i V_i^{-1} e_i$? No, the score needs $P e$ or directly $P V_k P y$. * Let $P_i e_i = V_i^{-1} e_i - V_i^{-1} X_i (X^T V^{-1} X)^{-1} X^T V^{-1} e_i$. * Ah! Let's compute $C = P y$. $P_i y_i = V_i^{-1} y_i - V_i^{-1} X_i \beta$. * Wait. $P y = V^{-1} y - V^{-1} X \hat{\beta}$. * So $C_i = V_i^{-1} e_i$. We already have this! * Score for $\phi_k$: $U_k = \frac{1}{2} \sum_i C_i^T V_{i,k} C_i - \frac{1}{2} \sum_i tr(P_i V_{i,k})$. * Wait, the score formula is: $U(\phi_k) = \frac{1}{2} y^T P V_k P y - \frac{1}{2} tr(P V_k)$. Since $C = P y$, $C_i = (P y)_i$. $y^T P V_k P y = \sum_i (P_i y_i)^T V_{i,k} (P_i y_i) = \sum_i C_i^T V_{i,k} C_i$. And $tr(P V_k) = \sum_i tr(P_i V_{i,k})$. * Great. 3. Compute $\mathcal{I}_{kl} = \frac{1}{2} \sum_i tr(P_i V_{i,k} P_i V_{i,l})$. 4. Compute step $\Delta \phi = \mathcal{I}^{-1} U$. 5. Apply step to $\sigma_u^2$, $\sigma_v^2$, $\phi$. 6. If $\phi$ goes outside (-1, 1), truncate. 7. If $\sigma$'s go negative, we can use log-parameterization or just truncate at 0. (Better to use log-parameterization in the FS algorithm, or simply truncate at 0 in the original scale and let FS handle it. Original `sae` truncates/adds small delta). Let's do FS on $\delta_u = \sigma_u^2$, $\delta_v = \sigma_v^2$, $\phi$. We can just do FS on the log scale to avoid boundary issues. * Let $\theta_1 = \log(\sigma_u^2)$, $\theta_2 = \log(\sigma_v^2)$, $\theta_3 = \text{atanh}(\phi)$. * Let $g = \text{diag}(\frac{\partial \phi_k}{\partial \theta_j})$. * The score for $\theta$ is $U_\theta = g \cdot U_\phi$. * The information for $\theta$ is $\mathcal{I}_\theta = g \cdot \mathcal{I}_\phi \cdot g + diag(U_\phi \cdot \frac{\partial^2 \phi}{\partial \theta^2})$ -> we can use the expected information ignoring the second term: $\mathcal{I}_\theta \approx g \cdot \mathcal{I}_\phi \cdot g$. * This is a standard Delta method transformation. We can simply perform FS on the transformed parameters. * Let's verify the gradient of the transformation. * $\frac{\partial \sigma_u^2}{\partial \theta_1} = \sigma_u^2$ -> $\rho$ (the autocorrelation) keep as $\rho$, or $\phi$. * Let's just use $\phi$ as a parameter with logit or tanh constraint. * FS in `sae` uses $\sigma_u^2$, $\sigma_v^2$, $\rho$. * Let's stick to the untransformed scale for $\sigma_u^2$, $\sigma_v^2$ and $\rho$ but truncate negative variance components and $|\rho| > 1$. If the step takes $\rho$ outside, it's truncated. * Actually, to be safe and fast, L-BFGS on the unconstrained transformed scale is easier to implement without the risk of Fisher Scoring diverging. Since we must compute the gradient anyway (or we can compute it exactly using Fisher scoring), we can feed it to L-BFGS, which is very robust and handles boundaries well via the transformation. * Let's backtrack. The question asks for an "efficient dan cepat" (efficient and fast) implementation. The FS algorithm is very fast if it converges. Implementing FS robustly is key. * Some small area estimation R packages with C++ backend use FS on the original parameters, with step halving if the likelihood decreases. This is called "Fisher Scoring with step control" (aka Newton-Raphson). * Let's do this. It's very standard. * Let's write the C++ code step by step. * **Class/Struct for Data:** * `int D; int T; int n; int p;` * `Eigen::MatrixXd X;` (n x p) * `Eigen::VectorXd y;` (n) * `Eigen::VectorXd vardir;` (n) * **Mappings:** * Index $i$ (domain) from $1$ to $D$, $t$ (time) from $1$ to $T$. * $k = (i-1)*T + (t-1)$. * **Function `LogLik`** (to check step qualities): * Calculate $\hat{\beta} = (X^T V^{-1} X)^{-1} X^T V^{-1} y$. * $l_{REML} = -\frac{1}{2} \log |V| - \frac{1}{2} \log |X^T V^{-1} X| - \frac{1}{2} (y - X\beta)^T V^{-1} (y - X\beta)$. * We can compute it to check if FS step increases likelihood. * **Function `FisherScoringStep`**: * Input: $\sigma_u^2, \sigma_v^2, \rho$. * Output: $\Delta \sigma_u^2, \Delta \sigma_v^2, \Delta \rho$. * Compute $V_i^{-1}$ using Cholesky. * Compute $\hat{\beta}$. * Compute $e_i = y_i - X_i \hat{\beta}$. * Compute $P_i e_i$ (which is $V_i^{-1} e_i - V_i^{-1} X_i (X^T V^{-1} X)^{-1} X^T V^{-1} e_i$). * Compute $U_{\sigma_u^2} = \frac{1}{2} \sum_i (P_i e_i)^T J (P_i e_i) - \frac{1}{2} \sum_i tr(P_i J)$. * Compute $U_{\sigma_v^2} = \frac{1}{2} \sum_i (P_i e_i)^T \Omega_i (P_i e_i) - \frac{1}{2} \sum_i tr(P_i \Omega_i)$. * Compute $U_{\rho} = \frac{1}{2} \sum_i (P_i e_i)^T \frac{\partial V_i}{\partial \rho} (P_i e_i) - \frac{1}{2} \sum_i tr(P_i \frac{\partial V_i}{\partial \rho})$. * Compute $\mathcal{I}_{\sigma_u^2, \sigma_u^2} = \frac{1}{2} \sum_i tr(P_i J P_i J)$. * Compute $\mathcal{I}_{\sigma_u^2, \sigma_v^2} = \frac{1}{2} \sum_i tr(P_i J P_i \Omega_i)$. * Compute $\mathcal{I}_{\sigma_u^2, \rho} = \frac{1}{2} \sum_i tr(P_i J P_i \frac{\partial V_i}{\partial \rho})$. * Compute $\mathcal{I}_{\sigma_v^2, \sigma_v^2} = \frac{1}{2} \sum_i tr(P_i \Omega_i P_i \Omega_i)$. * Compute $\mathcal{I}_{\sigma_v^2, \rho} = \frac{1}{2} \sum_i tr(P_i \Omega_i P_i \frac{\partial V_i}{\partial \rho})$. * Compute $\mathcal{I}_{\rho, \rho} = \frac{1}{2} \sum_i tr(P_i \frac{\partial V_i}{\partial \rho} P_i \frac{\partial V_i}{\partial \rho})$. * Solve $\mathcal{I} \Delta = U$. * Return $\Delta$. * **Function `compute_P_e`**: * $P y = V^{-1} y - V^{-1} X \hat{\beta}$. * $\hat{\beta} = (X^T V^{-1} X)^{-1} X^T V^{-1} y$. * Let $A = V^{-1} X$. ($n \times p$) * Let $B = X^T A$. ($p \times p$), the K matrix, already computed for $\hat{\beta}$. * $K = X^T V^{-1} X$. Cholesky: $R R^T = K$. * $C = V^{-1} y$. * $\hat{\beta} = K^{-1} X^T C$. * $P y = C - A \hat{\beta}$. * We need $P V_k P y$. * Wait, $P_i e_i = (C_i - A_i \hat{\beta})$. * Let's compute $w_i = P_i e_i = C_i - A_i \hat{\beta}$. * Score: $U_k = \frac{1}{2} \sum_i w_i^T V_{i,k} w_i - \frac{1}{2} \sum_i tr(P_i V_{i,k})$. * **Fisher Scoring loop:** * $\tau = 1$ (step size). * While $l_{REML}(\phi + \tau \Delta) < l_{REML}(\phi)$ and $\tau > 10^{-10}$: $\tau = \tau / 2$. * Update $\phi = \phi + \tau \Delta$. * Check convergence: $\max |\Delta| < \epsilon$ (e.g. $10^{-6}$). * **Bootstrap Function**: * Run FS on original data to get $\hat{\sigma}_u^2$, $\hat{\sigma}_v^2$, $\hat{\rho}$, $\hat{\beta}$. * Compute EBLUPs $\hat{\theta}_i$. * Pre-allocate results. * `#pragma omp parallel for` * For b in 1:B: * Generate $u_i \sim N(0, \hat{\sigma}_u^2)$ (vector size D). * Generate $v_{i1} \sim N(0, \hat{\sigma}_v^2 / (1-\hat{\rho}^2))$. * For $t=2$ to $T$: $v_{it} = \hat{\rho} v_{i,t-1} + \eta_{it}$, $\eta_{it} \sim N(0, \hat{\sigma}_v^2)$. * $\theta^b_i = X_i \hat{\beta} + u_i + v_i$. * $y^b_i = \theta^b_i + e_i$, $e_i \sim N(0, \psi_i)$. * Run FS on bootstrap data $(y^b, X, \psi)$ to get $\hat{\beta}^b$, $\hat{\sigma}_u^{2b}$, $\hat{\sigma}_v^{2b}$, $\hat{\rho}^b$. * Compute bootstrap EBLUP $\hat{\theta}^b_i$ using the FS estimated parameters. * Compute squared error $(\hat{\theta}^b_i - \theta^b_i)^2$. * Accumulate to `mse_sum`. * Return `mse = mse_sum / B`. * **Wait, let's double check the model** Does the model have a domain effect $u_i$ AND a domain-time AR1 effect $v_{it}$? Yes, the STFH model in `sae::pbmseSTFH`: $\theta_{it} = x_{it}^T \beta + u_i + v_{it}$ $v_{it} = \rho v_{i,t-1} + \epsilon_{it}$. $u_i \sim N(0, \sigma_u^2)$, $v_{i,1} \sim N(0, \sigma_v^2/(1-\rho^2))$, $\epsilon_{it} \sim N(0, \sigma_v^2)$. Direct estimator variance: $\psi_{it}$. * **Checking Efficiency Gains in C++** 1. **Matrix operations**: R's vectorized operations are already fast, but C++ avoids function call overhead and garbage collection, and allows OpenMP. 2. **OpenMP**: The bootstrap loop is embarrassingly parallel. We have `B` replicates. We can use `RcppParallel` or just standard OpenMP (`// [[Rcpp::plugins(openmp)]]`). 3. **Memory locality**: C++ can store and access data in contiguous vectors. 4. **Algorithmic optimizations**: * Pre-compute AR(1) matrices? The AR(1) matrix $\Omega_i(\rho)$ depends on $\rho$ and $T$. * We can precompute $\Omega_i$ basis matrices. $\Omega(\rho) = \frac{1}{1-\rho^2} C(\rho)$. * $C(\rho)$: $C_{ij} = \rho^{|i-j|}$. * We can compute $\frac{\partial V}{\partial \rho}$ during the Fisher scoring. * We can compute $V_i$ using a closed-form for the Cholesky of $V_i$ if it's purely AR1 + Diag(w). But the $u_i$ random effect adds $\sigma_u^2$ to every element of the block. This breaks pure AR1. * $V_i = \sigma_u^2 J_T + \frac{\sigma_v^2}{1-\rho^2} C + D_i$. * We can still compute $V_i$ as a $T \times T$ matrix and use Eigen's Cholesky. $T$ is usually small (e.g. $T <= 20$). $D$ can be very large (e.g. $D = 1000$). $D \times T$ can be 20000. * The REML loop loops over $i=1..D$. This is $O(D)$. * Inside the loop, $T$ is small. $T^3$ is tiny (8000 for T=20). * Bottleneck: the Fisher scoring loop is iterated over $D$ for every function evaluation. This is exactly what C++ excels at. R's loop overhead is huge here. * Profile the R function: `Rprof` would show the main time spent in the Fisher scoring iteration loop inside `REML` and the bootstrap loop. * Rewriting the Fisher scoring loop in C++ removes ALL R interpreter overhead from the inner loop. * **Let's write the C++ `REML` estimation function** We need a class or a set of functions. * `compute_V(mat, sig_u2, sig_v2, rho, vardir_i)`: * $\rho$? Let's call it `phi`. * $\Omega(\rho)$: known as `omega_mat(rho, T)`. * `V = sig_u2 * J + sig_v2 * omega + diag(vardir_i)`. * `compute_omega(rho, T)` * `Omega` matrix: $\frac{1}{1-\rho^2}$ `C_mat` where `C_mat(i,j) = rho^{|i-j|}`. * `compute_dV_drho(rho, T)`: * $\frac{\partial \Omega}{\partial \rho} = \frac{2\rho}{(1-\rho^2)^2} C + \frac{1}{1-\rho^2} C'$, where $C'(i,j) = |i-j| \rho^{|i-j|-1}$ (if $|i-j|>0$, else 0). * $\frac{\partial V}{\partial \rho} = \sigma_v^2 \frac{\partial \Omega}{\partial \rho}$. * `compute_PVPVtrace`: * $tr(P_i V_k P_i V_l) = tr( (P_i V_k) (P_i V_l) ) $. * $P_i = V_i^{-1} - V_i^{-1} X_i K^{-1} X_i^T V_i^{-1}$. * We can compute $C_i = V_i^{-1} y_i$. * $A_i = V_i^{-1} X_i$. (T x p matrix) * $K = \sum_i X_i^T A_i$. (p x p) * $R R^T = K$. * $\hat{\beta} = K^{-1} \sum_i X_i^T C_i$. * $w_i = C_i - A_i \hat{\beta}$. This is $P_i y_i$. * We need $tr(P_i V_{i,k})$. * $tr(P_i V_{i,k}) = tr(V_i^{-1} V_{i,k}) - tr(V_i^{-1} X_i K^{-1} X_i^T V_i^{-1} V_{i,k})$. * $tr(V_i^{-1} V_{i,k})$ is easy. * The second term: $V_{i,k}$ is symmetric. * $tr(K^{-1} X_i^T V_i^{-1} V_{i,k} V_i^{-1} X_i)$. * Let $B_{i,k} = X_i^T V_i^{-1} V_{i,k} V_i^{-1} X_i$. * Then the trace is $tr(K^{-1} \sum_i B_{i,k})$. * We need $tr(P_i V_k P_i V_l)$. * $P_i V_k = V_i^{-1} V_{i,k} - V_i^{-1} X_i K^{-1} X_i^T V_i^{-1} V_{i,k}$. * This is very messy to vectorize fully, but we can compute $P_i V_k$ as a $T \times T$ dense matrix and then $tr(P_i V_k P_i V_l)$ is just the Frobenius inner product of $P_i V_k$ and $(P_i V_l)^T$. Since $P_i V_l$ is symmetric, $tr(P_i V_k P_i V_l) = (P_i V_k) : (P_i V_l)$. We can compute this directly. * Wait, Rao & Molina (2015) Appendix C.2 has the exact formulas for Score and Fisher Information for the FH model with AR(1), using the spectral decomposition. * Actually, for $T \leq 10$, the full matrix multiplication in Eigen is already extremely fast. Let's just do the dense matrix multiplication. $O(D T^3)$ with $T=10$ is negligible. * Wait, we have a sum over $D$ domains. $D$ can be large. The $T \times T$ matrix operations are done $D$ times per Fisher scoring iteration. * Let's simplify the implementation of Score and Fisher Information. * Compute $V_i^{-1}$ (Cholesky). * Compute $w_i = V_i^{-1} e_i$ (T vector). * Compute $U_{\sigma_u^2}^i = \frac{1}{2} w_i^T J_T w_i - \frac{1}{2} tr(P_i J)$. * Compute $U_{\sigma_v^2}^i = \frac{1}{2} w_i^T \Omega_i w_i - \frac{1}{2} tr(P_i \Omega_i)$. * Compute $U_{\rho}^i = \frac{1}{2} w_i^T \frac{\partial V_i}{\partial \rho} w_i - \frac{1}{2} tr(P_i \frac{\partial V_i}{\partial \rho})$. * Trace of $P_i A$: $P_i = V_i^{-1} - V_i^{-1} X_i K^{-1} X_i^T V_i^{-1}$. $tr(P_i A) = tr(V_i^{-1} A) - tr(A V_i^{-1} X_i K^{-1} X_i^T V_i^{-1})$. $tr(V_i^{-1} A)$ is easy. $tr(A V_i^{-1} X_i K^{-1} X_i^T V_i^{-1}) = tr( K^{-1} X_i^T V_i^{-1} A V_i^{-1} X_i)$. Let $M_i = V_i^{-1} X_i$ (T x p). Let $N_i = V_i^{-1} A V_i^{-1} X_i$ (T x p). $tr(P_i A) = tr(V_i^{-1} A) - tr( K^{-1} X_i^T N_i )$. Since $K^{-1}$ is constant across domains, $\sum_i T_i = tr(K^{-1} \sum_i X_i^T N_i)$. * Let's apply this to $A = J$, $A = \Omega$, $A = \frac{\partial V}{\partial \rho}$. * Compute $N_i^J = V_i^{-1} J V_i^{-1} X_i$. * Compute $N_i^\Omega = V_i^{-1} \Omega V_i^{-1} X_i$. * Compute $N_i^{dV} = V_i^{-1} \frac{\partial V}{\partial \rho} V_i^{-1} X_i$. * Then $U_k = \frac{1}{2} \sum_i (w_i^T V_k w_i) - \frac{1}{2} ( \sum_i tr(V_i^{-1} V_k) - tr(K^{-1} \sum_i X_i^T N_i^k) )$. * For the Fisher Information $I_{kl} = \frac{1}{2} tr(P V_k P V_l)$: * $\sum_i tr(P_i V_k P_i V_l)$. * Let $B_i = V_k V_i^{-1} X_i$? No, $P_i V_k = V_i^{-1} V_k - V_i^{-1} X_i K^{-1} X_i^T V_i^{-1} V_k$. * This is $T \times T$ matrix multiplication. Let's just compute $P_i V_k$ and then trace against $P_i V_l$. * $P_i V_k = V_i^{-1} V_k - M_i K^{-1} M_i^T V_k$. * Compute $P_i V_k$ for all $k$ (3 matrices). Then $tr(P_i V_k P_i V_l)$. * Since $T$ is small, we can compute these as dense matrices. The overhead of writing out the complex trace formula vs. just doing dense matrix multiplication in Eigen is negligible. Let's just do dense matrix multiplication. It will be correct and fast. * **Simpler approach for Score and Info:** * `Eigen::MatrixXd Vi = V_i_inv;` (T x T) * `Eigen::MatrixXd Ji = J_T;` * `Eigen::MatrixXd Oi = Omega_i;` * `Eigen::MatrixXd dV = sig_v2 * dOmega;` * `Eigen::MatrixXd Pi = Vi - Vi * Xi * K_inv * Xi.transpose() * Vi;` * `Eigen::VectorXd wi = Vi * ei;` (Wait, $w_i = P_i y_i = V_i^{-1} e_i - M_i K^{-1} X_i^T V_i^{-1} e_i$. We can compute $w_i = Pi * y_i$). * Actually, let's just compute everything with matrix multiplication. It's $D$ times, $T$ is small. Running $D$ iterations of a Fisher scoring for $D=500, T=10$ in C++ is nearly instant. * $U_{\sigma^2_u} = \frac{1}{2} \sum_i w_i^T J_i w_i - \frac{1}{2} \sum_i tr(P_i J_i)$. * $tr(P_i J_i)$: `MatrixXd temp = Pi * Ji;` `double tr_Pi_Ji = temp.trace();` * $tr(P_i \Omega_i)$, $tr(P_i \frac{\partial V_i}{\partial \rho})$ are similar. * $I_{kl} = \frac{1}{2} \sum_i tr(P_i V_k P_i V_l)$. `MatrixXd temp1 = Pi * Ji;` `MatrixXd temp2 = Pi * Oi;` `double I_sig_sig = 0.5 * (temp1 * temp2.transpose()).trace();` Wait, $tr(A B) = (A^T) : (B)$. If we compute $C = P_i J$, $D = P_i \Omega$. $tr(C D) = (C^T : D) = vec(C)^T vec(D)$. In Eigen, `(Pi*Ji).cwiseProduct((Pi*Oi).transpose()).sum()` or just `(Pi*Ji * Pi*Oi).trace()`. Yes, $tr(M_1 M_2) = \sum_{i} \sum_{j} (M_1)_{ij} (M_2)_{ji}$. `(M1 * M2).trace()` computes this. So $I_{kl} = 0.5 * \sum_i ( (P_i V_k) * (P_i V_l) ).trace()$. * This is very clean to code. We just need to manage the sums over domains. * We can do this in a single loop over domains. * **Data Storage** * `y_i`: a vector of length $T$ for each domain $i$. We can store the entire `y` as an `Eigen::VectorXd` of length $n = D * T$, or an `Eigen::MatrixXd` of size $T \times D$. * Let's store it as a vector for simplicity. * `X_i`: matrix $T \times p$. * `vardir_i`: vector of length $T$. * To extract: `Eigen::Map yi(y.data() + i * T, T);` `Eigen::MatrixXd Xi = X.block(p * i, 0, T, p);` Wait, if `X` is a matrix, we can do `X.middleRows(i * T, T)`. * **Optimizer Details** * We will write a simple Fisher Scoring function with step halving. * `int reml_fs(Eigen::VectorXd& params, const Eigen::VectorXd& y, const Eigen::MatrixXd& X, const Eigen::VectorXd& vardir, int D, int T, int max_iter = 100, double tol = 1e-6)` * `params` is a vector of size 3: `(sig_u2, sig_v2, rho)`. * Return value: 0 for success, -1 for failure. * Inside, we compute `V_i_i`, `V_i_omega_i`, `V_i_dV_id_rho`, etc. We will need to compute $\Omega_i(\rho)$ for a given $\rho$. * $\Omega_i$ matrix: $\Omega_i(\rho) = \frac{1}{1-\rho^2} \begin{pmatrix} 1 & \rho & \rho^2 & \dots \\ \rho & 1 & \rho & \dots \\ \vdots & \vdots & \ddots \end{pmatrix}$ * $\frac{\partial \Omega_i}{\partial \rho} = \frac{2\rho}{(1-\rho^2)^2} \Omega_i + \frac{1}{1-\rho^2} \frac{\partial C}{\partial \rho}$. Wait, standard AR(1) variance matrix: $\text{Var}(v) = \frac{\sigma_v^2}{1-\rho^2} \begin{pmatrix} 1 & \rho & \dots \\ \rho & 1 & \dots \\ \vdots & \vdots & \ddots \end{pmatrix}$. Let's just compute $\Sigma_v = \frac{\sigma_v^2}{1-\rho^2} C$. $C_{st} = \rho^{|s-t|}$. This is easy to construct: ```cpp Eigen::MatrixXd compute_C(double rho, int T) { Eigen::MatrixXd C(T, T); for (int i = 0; i < T; ++i) { for (int j = 0; j < T; ++j) { C(i,j) = std::pow(rho, std::abs(i - j)); } } return C; } Eigen::MatrixXd compute_V() { return sig_u2 * Eigen::MatrixXd::Ones(T, T) + sig_v2 / (1.0 - rho*rho) * C + vardir_i.asDiagonal().toDenseMatrix(); } ``` * Actually, the standard reparameterization to make it scale free: Let $\sigma_v^2$ be the variance of the transition errors. Var(v_1) = $\sigma_v^2/(1-\rho^2)$. Var(v) = $\sigma_v^2/(1-\rho^2) * C$. * **Corner Cases / Convergence Issues:** 1. $\rho$ close to 1 or -1. 2. $\sigma_u^2$ or $\sigma_v^2$ close to 0. 3. Fisher scoring overshoots. * Step halving: `while (new_lml > old_lml + step_tol)` halve step. * If $|rho| > 0.995$, truncate to 0.995. * If sigma's go negative, set to a small positive value (e.g. `1e-6`). * **Parallel Bootstrap Loop** ```cpp #include ... struct BootstrapWorker : public RcppParallel::Worker { ... void operator()(std::size_t begin, std::size_t end) { for (std::size_t b = begin; b < end; ++b) { // Generate theta_b, y_b // Run REML // Compute EBLUP // Accumulate } } }; ``` Or simpler, just use `#pragma omp parallel for` if available. Since many R packages do this, it's safe. ```cpp #pragma omp parallel for for (int b = 0; b < B; ++b) { ... } ``` * **Wait, the R package `sae` already uses C?** Let's check the source code of `sae` on R-Forge or GitHub. The `pbmseSTFH` is an R wrapper that calls internal functions. "The STFH model is fitted by Fisher-scoring algorithm and bootstrap mean squared error is obtained." The internal functions `stfh` are written in R with some calls to Fortran/C for matrix operations? Let's check the source: `sae` package. The main functions for `pbmseSTFH` are in R. They call `stFH` or similar which is an R function. The R function `stFH` does the Fisher scoring loop in R! This means the bootstrap loop in R calls R function `stFH`. Wait, if the R function does the FS loop in R, it's extremely slow because R loops are slow. *Translating the FS loop to C++ makes it extremely fast.* The difference between R's FS loop and C++'s FS loop is a classic benchmark (100x speedup on the loop itself is common). * **Let's explicitly fix the model to exactly match `sae::pbmseSTFH`.** * Model: $\theta_{it} = x_{it}^T \beta + u_i + v_{it}$. * $u_i \sim N(0, \sigma_u^2)$. * $v_{it} = \rho v_{i,t-1} + \epsilon_{it}$, $\epsilon_{it} \sim N(0, \sigma_v^2)$. * $v_{i,1} \sim N(0, \sigma_v^2/(1-\rho^2))$. * Estimator: EBLUP. $\hat{\theta}_{it} = x_{it}^T \hat{\beta} + \hat{u}_i + \hat{v}_{it}$. * MSE: Parametric bootstrap. * **Generate Bootstrap Sample:** 1. Get estimated parameters: $\hat{\beta}_{GLS}, \hat{\sigma}_u^2, \hat{\sigma}_v^2, \hat{\rho}$ from the original sample. 2. Calculate $\hat{\theta}$ (original EBLUP). 3. For b in 1:B: * Draw $u_i^{(b)} \sim N(0, \hat{\sigma}_u^2)$. * Draw $v_{i,1}^{(b)} \sim N(0, \hat{\sigma}_v^2/(1-\hat{\rho}^2))$. * Draw $\epsilon_{it}^{(b)} \sim N(0, \hat{\sigma}_v^2)$. * $v_{it}^{(b)} = \hat{\rho} v_{i,t-1}^{(b)} + \epsilon_{it}^{(b)}$. * $\theta_{it}^{(b)} = x_{it}^T \hat{\beta} + u_i^{(b)} + v_{it}^{(b)}$. * Draw $e_{it}^{(b)} \sim N(0, \psi_{it})$ (sampling error). * $y_{it}^{(b)} = \theta_{it}^{(b)} + e_{it}^{(b)}$. * Fit the model to $y^{(b)}$: get $\hat{\beta}^{(b)}, \hat{\sigma}_u^{2(b)}, \hat{\sigma}_v^{2(b)}, \hat{\rho}^{(b)}$. * Compute $\hat{\theta}_{it}^{(b)}$ (bootstrap EBLUP). * Loop. * **Could we reuse the R random number generator?** Yes, we should use R's RNG for reproducibility. ```cpp #include Rcpp::RNGScope scope; // Not strictly needed in modern Rcpp // Use R::rnorm() or Rcpp::rnorm() ``` * **Detailed C++ Pseudo-code** ```cpp #include #include // [[Rcpp::depends(RcppEigen)]] // [[Rcpp::depends(RcppParallel)]] // [[Rcpp::plugins(openmp)]] using namespace Rcpp; using Eigen::MatrixXd; using Eigen::VectorXd; using Eigen::LLT; // Compute auto-correlation matrix C: C_{ij} = rho^{|i-j|} MatrixXd ar1_corr(double rho, int T) { MatrixXd C(T, T); for (int i = 0; i < T; ++i) for (int j = 0; j < T; ++j) C(i,j) = std::pow(rho, std::abs(i - j)); return C; } // Compute derivative of C w.r.t rho: dC_{ij} = |i-j| * rho^{|i-j|-1} MatrixXd ar1_dcorr(double rho, int T) { MatrixXd dC(T, T); double eps = (rho == 0 && std::abs(rho) < 1e-10) ? 0 : 1; // Handle limit when rho -> 0 for (int i = 0; i < T; ++i) { for (int j = 0; j < T; ++j) { int lag = std::abs(i - j); if (lag == 0) { dC(i,j) = 0.0; } else if (lag == 1 && rho == 0) { dC(i,j) = 1.0; // derivative of rho^1 = 1 } else if (rho == 0) { dC(i,j) = 0.0; } else { dC(i,j) = lag * std::pow(rho, lag - 1); } } } return dC; } // Compute -2 * Restricted Log-Likelihood double neg2_log_lik(const VectorXd& params, const VectorXd& y, const MatrixXd& X, const VectorXd& vardir, int D, int T) { double sig_u2 = params[0]; double sig_v2 = params[1]; double rho = params[2]; if (sig_u2 <= 0 || sig_v2 <= 0 || std::abs(rho) >= 1.0) return 1e15; if (sig_v2 <= 0) return 1e15; double val = 0.0; int n = D * T; int p = X.cols(); // Compute J_T (matrix of ones) MatrixXd J_T = MatrixXd::Constant(T, T, 1.0); MatrixXd C = ar1_corr(rho, T); MatrixXd InvSigmaSum = MatrixXd::Zero(p, p); VectorXd XtSigma_y = VectorXd::Zero(p); double log_det_V = 0.0; bool first = true; for (int d = 0; d < D; ++d) { VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); // Build V_d MatrixXd V = sig_u2 * J_T + sig_v2 / (1.0 - rho*rho) * C + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); if (lltV.info() != Eigen::Success) return 1e15; MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); double log_det_V_d = 2.0 * lltV.matrixL().toDenseMatrix().diagonal().array().log().sum(); log_det_V += log_det_V_d; // X^T V^{-1} X MatrixXd XtVinv = Xd.transpose() * V_inv; InvSigmaSum += XtVinv * Xd; XtSigma_y += XtVinv * yd; } // Beta_hat LLT lltX(InvSigmaSum); if (lltX.info() != Eigen::Success) return 1e15; VectorXd beta = lltX.solve(XtSigma_y); double log_det_XtVX = 2.0 * lltX.matrixL().toDenseMatrix().diagonal().array().log().sum(); double res_quad = 0.0; for (int d = 0; d < D; ++d) { VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2 * J_T + sig_v2 / (1.0 - rho*rho) * C + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); VectorXd resid = yd - Xd * beta; VectorXd Vinv_resid = V_inv * resid; res_quad += resid.dot(Vinv_resid); } val = log_det_V + log_det_XtVX + res_quad; return val; } // Main Fisher Scoring function bool fit_stfh_reml(VectorXd& params, const VectorXd& y, const MatrixXd& X, const VectorXd& vardir, int D, int T, int max_iter=100, double tol=1e-6) { double sig_u2 = params[0]; double sig_v2 = params[1]; double rho = params[2]; int p = X.cols(); MatrixXd J_T = MatrixXd::Constant(T, T, 1.0); MatrixXd Eye = MatrixXd::Identity(T, T); for (int iter = 0; iter < max_iter; ++iter) { MatrixXd C = ar1_corr(rho, T); MatrixXd dC = ar1_dcorr(rho, T); double inv_den = 1.0 / (1.0 - rho*rho); double d_inv_den = 2.0 * rho / ((1.0 - rho*rho) * (1.0 - rho*rho)); MatrixXd Omega = inv_den * C; MatrixXd dOmega = d_inv_den * C + inv_den * dC; VectorXd U = VectorXd::Zero(3); MatrixXd I = MatrixXd::Zero(3, 3); MatrixXd InvSigmaSum = MatrixXd::Zero(p, p); VectorXd XtSigma_y = VectorXd::Zero(p); // First pass: compute beta_hat for (int d = 0; d < D; ++d) { VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2 * J_T + sig_v2 * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(Eye); MatrixXd XtVinv = Xd.transpose() * V_inv; InvSigmaSum += XtVinv * Xd; XtSigma_y += XtVinv * yd; } LLT lltK(InvSigmaSum); if (lltK.info() != Eigen::Success) { // params = ... // try steepest descent or something // Rcout << "lltK failed" << std::endl; return false; } VectorXd beta = lltK.solve(XtSigma_y); MatrixXd K_inv = lltK.solve(Eye); // Second pass: compute Score and Info for (int d = 0; d < D; ++d) { VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2 * J_T + sig_v2 * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(Eye); VectorXd resid = yd - Xd * beta; VectorXd Vinv_resid = V_inv * resid; MatrixXd XtVinv = Xd.transpose() * V_inv; // P = V^{-1} - V^{-1} X K^{-1} X^T V^{-1} MatrixXd P = V_inv - V_inv * Xd * K_inv * XtVinv; // V_k: k=0 sig_u2 -> J_T // k=1 sig_v2 -> Omega // k=2 rho -> sig_v2 * dOmega MatrixXd V_k[3] = {J_T, Omega, sig_v2 * dOmega}; for (int k = 0; k < 3; ++k) { VectorXd w_k = V_k[k] * Vinv_resid; U(k) += 0.5 * (Vinv_resid.dot(w_k) - (P * V_k[k]).trace()); for (int l = k; l < 3; ++l) { MatrixXd PV_k = P * V_k[k]; MatrixXd PV_l = P * V_k[l]; I(k,l) += 0.5 * (PV_k * PV_l).trace(); } } } // Symmetrize I for (int k = 0; k < 3; ++k) for (int l = 0; l < k; ++l) I(l,k) = I(k,l); // Solve I * delta = U LLT lltI(I); if (lltI.info() != Eigen::Success) { // Maybe a small diagonal perturbation for stability I += 1e-8 * MatrixXd::Identity(3,3); lltI.compute(I); if (lltI.info() != Eigen::Success) return false; } VectorXd delta = lltI.solve(U); // Step halving double step_size = 1.0; VectorXd params_new = params; double old_neg2ll = neg2_log_lik(params, y, X, vardir, D, T); for (int h = 0; h < 20; ++h) { params_new = params + step_size * delta; // Ensure positivity and |rho| < 1 if (params_new[0] < 0) params_new[0] = 1e-8; if (params_new[1] < 0) params_new[1] = 1e-8; if (std::abs(params_new[2]) >= 1.0) { if (params_new[2] >= 1.0) params_new[2] = 0.99; if (params_new[2] <= -1.0) params_new[2] = -0.99; } double new_neg2ll = neg2_log_lik(params_new, y, X, vardir, D, T); if (std::isnan(new_neg2ll) || new_neg2ll > old_neg2ll + 1e-6) { step_size *= 0.5; } else { break; } } params = params_new; if (delta.norm() < tol) break; } return true; } // [[Rcpp::export]] Rcpp::List pbmseSTFH_fast(Rcpp::NumericVector yR, Rcpp::NumericMatrix XR, int D, int T, Rcpp::NumericVector vardirR, int B, int max_iter = 100) { int n = D * T; int p = XR.ncol(); // Map to Eigen Eigen::Map y(yR.begin(), n); Eigen::Map X(XR.begin(), n, p); Eigen::Map vardir(vardirR.begin(), n); // Initial values VectorXd params(3); params[0] = 0.5 * vardir.mean(); // sig_u2 params[1] = 0.5 * vardir.mean(); // sig_v2 params[2] = 0.5; // rho // Fit original model fit_stfh_reml(params, y, X, vardir, D, T, max_iter); double sig_u2_hat = params[0]; double sig_v2_hat = params[1]; double rho_hat = params[2]; // Compute original EBLUP VectorXd theta_hat_orig(n); { MatrixXd J_T = MatrixXd::Constant(T, T, 1.0); MatrixXd C = ar1_corr(rho_hat, T); double inv_den = 1.0 / (1.0 - rho_hat*rho_hat); MatrixXd Omega = inv_den * C; MatrixXd InvSigmaSum = MatrixXd::Zero(p, p); VectorXd XtSigma_y = VectorXd::Zero(p); for (int d = 0; d < D; ++d) { VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2_hat * J_T + sig_v2_hat * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); MatrixXd XtVinv = Xd.transpose() * V_inv; InvSigmaSum += XtVinv * Xd; XtSigma_y += XtVinv * yd; } LLT lltK(InvSigmaSum); VectorXd beta_hat = lltK.solve(XtSigma_y); MatrixXd K_inv = lltK.solve(MatrixXd::Identity(p, p)); for (int d = 0; d < D; ++d) { VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2_hat * J_T + sig_v2_hat * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); MatrixXd XtVinv = Xd.transpose() * V_inv; // Gamma = Z^T V^{-1} Z // Z = [1_T | I_T] // Cov matrix of (u_i, v_i) = diag(sig_u2, sig_v2 * Omega) // Actually the EBLUP formula: \hat{\theta}_i = X_i \hat{\beta} + G_i Z_i^T V_i^{-1} (y_i - X_i \hat{\beta}) // Z_i = [1_T, I_T] // G_i = blkdiag(sig_u2, sig_v2 * Omega) // Z_i^T = [1_T^T ; I_T] // G_i Z_i^T = [sig_u2 * 1_T^T ; sig_v2 * Omega] // Let's compute GZ = [sig_u2 * 1_T, sig_v2 * Omega] (T x (1+T)) VectorXd resid = yd - Xd * beta_hat; VectorXd Vinv_resid = V_inv * resid; // u_i_hat = sig_u2 * 1_T^T * Vinv_resid double u_hat = sig_u2_hat * VectorXd::Ones(T).dot(Vinv_resid); // v_i_hat = sig_v2_hat * Omega * Vinv_resid VectorXd v_hat = sig_v2_hat * Omega * Vinv_resid; // In the STFH model, u_i is domain effect, v_i is domain-time effect. // EBLUP: \hat{\theta}_i = X_i \hat{\beta} + u_i_hat + v_i_hat VectorXd Xb = Xd * beta_hat; theta_hat_orig.segment(d * T, T) = Xb + u_hat * VectorXd::Ones(T) + v_hat; } } // Bootstrap VectorXd mse_sum = VectorXd::Zero(n); #pragma omp parallel for for (int b = 0; b < B; ++b) { VectorXd y_b(n); VectorXd theta_b(n); // Generate bootstrap data for (int d = 0; d < D; ++d) { VectorXd psi = vardir.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); double u_b = R::rnorm(0, std::sqrt(sig_u2_hat)); VectorXd v_b(T); // Generate AR(1) process v_b[0] = R::rnorm(0, std::sqrt(sig_v2_hat / (1.0 - rho_hat*rho_hat))); for (int t = 1; t < T; ++t) { v_b[t] = rho_hat * v_b[t-1] + R::rnorm(0, std::sqrt(sig_v2_hat)); } VectorXd theta_linear = Xd * beta_hat; theta_b.segment(d * T, T) = theta_linear + u_b + v_b; for (int t = 0; t < T; ++t) { double e_b = R::rnorm(0, std::sqrt(psi[t])); y_b[d * T + t] = theta_b[d * T + t] + e_b; } } // Fit REML on bootstrap data VectorXd params_b(3); params_b[0] = sig_u2_hat; params_b[1] = sig_v2_hat; params_b[2] = rho_hat; fit_stfh_reml(params_b, y_b, X, vardir, D, T, 50); // less iterations // Calculate bootstrap EBLUP VectorXd theta_b_hat(n); { double sig_u2_b = params_b[0]; double sig_v2_b = params_b[1]; double rho_b = params_b[2]; MatrixXd C = ar1_corr(rho_b, T); double inv_den = 1.0 / (1.0 - rho_b * rho_b); MatrixXd Omega = inv_den * C; MatrixXd InvSigmaSum_b = MatrixXd::Zero(p, p); VectorXd XtSigma_y_b = VectorXd::Zero(p); for (int d = 0; d < D; ++d) { VectorXd yd_b = y_b.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2_b * MatrixXd::Constant(T, T, 1.0) + sig_v2_b * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); MatrixXd XtVinv = Xd.transpose() * V_inv; InvSigmaSum_b += XtVinv * Xd; XtSigma_y_b += XtVinv * yd_b; } LLT lltK_b(InvSigmaSum_b); VectorXd beta_b = lltK_b.solve(XtSigma_y_b); for (int d = 0; d < D; ++d) { VectorXd yd_b = y_b.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2_b * MatrixXd::Constant(T, T, 1.0) + sig_v2_b * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); // Compute EBLUP exactly as before for the bootstrap VectorXd resid_b = yd_b - Xd * beta_b; VectorXd Vinv_resid_b = V_inv * resid_b; double u_b_hat = sig_u2_b * VectorXd::Ones(T).dot(Vinv_resid_b); VectorXd v_b_hat = sig_v2_b * Omega * Vinv_resid_b; VectorXd Xb = Xd * beta_b; theta_b_hat.segment(d * T, T) = Xb + u_b_hat * VectorXd::Ones(T) + v_b_hat; } } VectorXd diff = theta_b_hat - theta_b; for (int i = 0; i < n; ++i) { #pragma omp atomic mse_sum[i] += diff[i] * diff[i]; } } VectorXd mse = mse_sum / B; return Rcpp::List::create(Rcpp::Named("mse") = Rcpp::wrap(mse), Rcpp::Named("est") = Rcpp::wrap(theta_hat_orig)); } ``` * **Wait, there might be a bug in `fit_stfh_reml` regarding the score and information calculation in the same loop vs separate loops for `beta`.** Let's double check the formula for $U(\rho)$ and $I(\rho, \rho)$. $U(\rho) = \frac{1}{2} y^T P \dot{V} P y - \frac{1}{2} tr(P \dot{V})$. We computed `w_k = V_k[2] * Vinv_resid`. Wait. `VectorXd w_k = V_k[k] * Vinv_resid;` $y^T P V_k P y = (P y)^T V_k (P y) = w^T V_k w$. Since we computed $w = P y = Vinv_resid - V_inv Xd K^{-1} XtVinv y$. In my loop, I computed $w_k = V_k * (y^T P)$? No. $P y = V^{-1} (y - X \beta) = Vinv_resid$. Wait! $P y = V^{-1} y - V^{-1} X \hat{\beta}$. $V^{-1} y - V^{-1} X \beta = V^{-1} (y - X\beta) = V^{-1} e$. So $P y = V^{-1} e$ ONLY if $\beta$ is the GLS estimator. Yes, $\hat{\beta} = (X^T V^{-1} X)^{-1} X^T V^{-1} y$. $P y = V^{-1} y - V^{-1} X \hat{\beta} = V^{-1} (y - X \hat{\beta}) = V^{-1} e$. So $w = P y = V^{-1} e$. Therefore $w^T V_k w = e^T V^{-1} V_k V^{-1} e$. This simplifies the score code! No, wait. We computed `w_k = V_k[k] * Vinv_resid`. $y^T P V_k P y = w^T V_k w$. In my loop I did `VectorXd w_k = V_k[k] * Vinv_resid;` `U(k) += 0.5 * (Vinv_resid.dot(w_k) - (P * V_k[k]).trace());` This computes $w^T V_k w = (Vinv\_resid)^T V_k (Vinv\_resid)$. Yes, `Vinv_resid.dot(w_k)` is exactly $w^T V_k w$. Trace of $P V_k$: `(P * V_k[k]).trace()` is $tr(P V_k)$. So the formula is correct! Let's check $I_{kl} = \frac{1}{2} tr(P V_k P V_l)$. `MatrixXd PV_k = P * V_k[k];` `MatrixXd PV_l = P * V_k[l];` `I(k, l) += 0.5 * (PV_k * PV_l).trace();` $tr(P V_k P V_l) = tr( (P V_k) (P V_l) ) = \sum_i \sum_j (P V_k)_{ji} (P V_l)_{ij}$. The trace of `PV_k * PV_l` is exactly $tr(P V_k P V_l)$. The formula is correct! *Important nuance*: The Score and Info depend on $\beta$ which depends on the current $V$. We computed a single $\beta$ for the current $V$, then looped to compute Score and Info. This is exactly correct for the REML Fisher scoring algorithm. * **Check the Step Halving function `neg2_log_lik`** The step halving requires computing the likelihood. I already wrote `neg2_log_lik`. It's computationally heavy ($O(D T^3)$). It's called ~20 times if the step is bad. Since $T$ is small, this is acceptable. If $T$ is large, we might want an embedded step-halving check using the Score ($U^T \Delta < 0$ or similar). But the likelihood check is the gold standard. * **Memory and OpenMP** Each thread needs its own copy of $y_b, X, vardir, params$. $n$ can be large. Making copies of $X$ for each thread is a problem. We should avoid copies of $X$ inside the OpenMP loop. Instead, we pass $X$ as a const reference. OpenMP sharing is fine for read-only data. ```cpp const Eigen::MatrixXd& X_ref = X; const Eigen::VectorXd& vardir_ref = vardir; const Eigen::VectorXd& y_const = y; // not used inside loop directly, y_b is used ``` Wait, `fit_stfh_reml` takes `const MatrixXd& X, const VectorXd& vardir`. These will be shared across threads. Good. But `fit_stfh_reml` modifies the data? No, it doesn't. It reads `y_b`. The OpenMP loop correctly operates on shared `X`, `vardir`, `beta_hat`, `sig_u2_hat`, etc., and private `y_b`, `params_b`. * **Wait, `#pragma omp atomic` might be very slow for a vector `mse_sum`.** `for (int i = 0; i < n; ++i) { #pragma omp atomic mse_sum[i] += ... }` This is a sequential bottleneck. If `n` is large, atomic operations on each element of `mse_sum` can be slow. *Alternative:* Use a reduction on the entire vector? OpenMP <= 4.0 does not support user-defined reductions of custom objects easily in all compilers. *Alternative 2:* Use a local `VectorXd mse_local` for each thread, initialized to 0. Accumulate into `mse_local`. At the end of the parallel block, combine into `mse_sum` using critical or atomic. Actually, each thread can have its own `mse_local`. ```cpp std::vector mse_local_vec(omp_get_max_threads(), VectorXd::Zero(n)); #pragma omp parallel for for (int b = 0; b < B; ++b) { int tid = omp_get_thread_num(); VectorXd diff = ...; mse_local_vec[tid] += diff.cwiseAbs2(); // or cwiseProduct } for (int t = 0; t < omp_get_max_threads(); ++t) { mse_sum += mse_local_vec[t]; } ``` *Wait!* `mse_local_vec[tid] += diff.cwiseAbs2()` is not thread safe if multiple threads access the same element of the map! Yes, it is. Each thread has its own `VectorXd mse_local_vec`. They don't overlap. Let's fix the memory for `mse_local_vec`. This is highly optimized. * **RNG in OpenMP** `R::rnorm()` is not thread safe! We must use a thread-safe RNG, or lock around `R::rnorm` (which defeats the purpose). ``R::rnorm()` calls internal R C API which is not thread safe. * *Solution:* Use C++11 `` in each thread, or the RcppParallel RNG wrapper, or `Rcpp::rnorm()` wrapper that might be thread-safe? No, `Rcpp::rnorm` is not thread safe either. *Best solution for R and Rcpp:* Use `Rcpp::rnorm()` outside the parallel loop to generate all bootstrap errors, or use a thread-safe RNG. Since the entire point is speed, we can generate all the bootstrap data and store it, then do the bootstrap loop in parallel? No, memory is $B \times n$ which could be huge. *Alternative:* Use `dqrng` or `Rcpp`'s `RNGScope`? *Actually, `R::rnorm()` can be called inside an OpenMP parallel region, but the RNG state is not protected. It can cause crashes or duplicates.* *Better:* Use `Rcpp::rnorm()` and the RNG is locked. Wait, the RNG scope can be serialized. The simplest and safest approach for high performance in R is to use `RcppZiggurat` or just use `R::runif`/`R::rnorm` from the thread-safe `Rmath` library? No, `R::rnorm` calls the internal R RNG which is NOT thread safe. Let's use `RcppParallel` RNG or just `std::mt19937`. *Alternative: Pre-generate the random numbers in R, pass them to C++.* This is the safest and easiest to code. ```R # In R: boot_u <- matrix(rnorm(D * B, 0, sqrt(sig_u2_hat)), ncol = B) ... ``` Passing a `Rcpp::NumericMatrix` of random draws is very fast and avoids all thread-safety issues. Let's design the function to accept pre-generated random draws. Wait, the bootstrap sample depends on `u`, `v`, `e`. The EBLUP depends on `y_b`. Passing all of them solves the RNG issue. Let's pass matrices: `Rcpp::NumericMatrix u_boot(D, B);` `Rcpp::NumericMatrix v_boot_wil(T, B)?` No, v_boot depends on domain and time. $e_{it}$ is $n \times B$. Actually, generating the random numbers in the C++ code using a parallelizable RNG like `dqrng` (which is thread safe) is better. Let's just use `Rcpp::rnorm()` in a single threaded loop BEFORE the parallel loop? This might be sequential. Since $B$ can be 500, generating $n * B$ random numbers sequentially is a very small overhead compared to the Fisher Scoring. Let's generate all random numbers in a single sequential loop, then run the bootstrap in parallel. Structure: 1. Generate `u_errors` (D x B matrix), `v_errors` (D x T x B array? but in C++, we can generate a flat vector of length `D * T * B` for the AR1 innovations, and `D * B` for $u_i$ and $v_{i1}$). 2. Put them in an Eigen matrix: `Eigen::MatrixXd u_b(D, B);` `Eigen::MatrixXd e_b(n, B);` `Eigen::MatrixXd eta_b(D, T-1, B)`? Let's keep it simple and generate them on the fly inside the OpenMP loop using a thread-safe RNG. `std::mt19937 rng(omp_get_thread_num() * time(0));` (or `std::random_device`). This avoids the memory overhead of storing B copies of everything, which might be 500 * 20000 = 10 million doubles, okay maybe it's fine, but it scales with D and T. Let's avoid storing them. Thread safe RNG in OpenMP: ```cpp #include #pragma omp parallel { std::mt19937 rng(omp_get_thread_num() + std::random_device{}()); std::normal_distribution norm(0, 1); #pragma omp for for (int b = 0; b < B; ++b) { // generate } } ``` This is purely C++, it doesn't touch R's RNG. It is very fast. * **Let's refine the `fit_stfh_reml` function.** Is there any speed bottleneck? The `neg2_log_lik` function is called in step halving. It loops over all domains. It repeats many computations. We can compute the log-lik more efficiently inside the main Fisher Scoring loop. But since step halving is rare (when the step is good, step=1), the 20 iterations are rarely fully executed. If the current $\beta$ is close, the step is usually good. So the overhead of computing the log-lik is acceptable. One more thing: `LLT` decomposition of $V$ each time. $V$ is $T \times T$. $T$ max maybe 20. This is extremely fast. $D$ up to 5000. So $D * 20^3$ is $5000 * 8000 = 40,000,000$ operations. This is trivial in C++. Let's profile the mental execution. The loop inside `fit_stfh_reml` calls: `ar1_corr(rho)` (T x T, O(T^2)) `ar1_dcorr(rho)` (T x T, O(T^2)) Then two loops over D: 1st loop: Build V, Cholesky, V^-1, X^T V^-1 X. (O(D T^3)) 2nd loop: Build V, Cholesky, V^-1, P, V_k, trace. (O(D k^2 T^3), k=3). This is extremely efficient. The whole REML for 500 domains should take microseconds. Bootstrap: B times the above. **What if $D$ is small but $T$ is large?** $T$ could be 100. $T^3$ is 1e6. $D \times T^3$ is $500 \times 1e6 = 5e8$. This is heavy for 1 iteration, let alone 100 iterations * B. If $T$ is large, the standard AR(1) Cholesky is not fully efficient (banded matrix). However, the STFH model adds $\sigma_u^2 J$, which makes it dense. If $\sigma_u^2$ is zero, $V$ is just an AR(1) + diagonal. The inverse of AR(1) + diagonal has a known $O(T)$ Kalman Filter form. If $\sigma_u^2$ is non-zero, $V$ is dense. We can still use the Kalman Filter for the full model? Yes, the STFH model ($u_i + v_{it}$) can be written as a state-space model and computed in $O(T)$ using the Kalman filter. This is the standard state-space approach (e.g. `STSAE`, or Rao & Molina appendix). *However*, for the typical use case, $T$ is small to moderate (e.g., 10-20 years of data). For $T=100$, the state-space approach would be much faster. Let's provide a *general* solution using the $T \times T$ matrix approach, and mention that if $T$ is large, a Kalman filter formulation is advisable. (The user can adapt it). * **Let's check for any missing elements in the standard algorithm** Standard `pbmseSTFH` returns `mse` and `est` (EBLUP)? Yes, and also `components` (random effects)? Let's just return `mse` and `est`. * **Wait, the bootstrap REML estimation might fail for some generated samples.** We should ensure robustness. If `fit_stfh_reml` returns `false`, we can skip the contribution for that bootstrap replicate, or set the mse to something, or use the previous bootstrap parameters. What does the original `pbmseSTFH` do? It probably has a `tryCatch`. If REML fails, the bootstrap replicate is discarded? In R, it is discarded and a warning is printed. In C++, we can count how many failed and adjust the MSE denominator. ```cpp int n_failed = 0; #pragma omp parallel for reduction(+:n_failed) for (int b = 0; b < B; ++b) { bool success = fit_stfh_reml(...); if (!success) { n_failed++; continue; } // calculate EBLUP } // mse = mse_sum / (B - n_failed); ``` * **Let's integrate the RNG and OpenMP properly.** ```cpp #include #include #include #ifdef _OPENMP #include #endif // [[Rcpp::depends(RcppEigen)]] // [[Rcpp::plugins(openmp)]] // ... functions ... // [[Rcpp::export]] Rcpp::List pbmseSTFH_fast(...) { // ... setup ... int n_failed = 0; #pragma omp parallel reduction(+:n_failed) { std::mt19937 rng(omp_get_thread_num() + 1); std::normal_distribution norm(0, 1); #pragma omp for for (int b = 0; b < B; ++b) { // generate bootstrap data using norm(rng) // fit model // compute EBLUP // update mse_local } } // combine mse_local // compute mse ... } ``` Wait, `std::mt19937` and `std::normal_distribution` are inside the parallel region. The `rng` is local to each thread. This is the standard pattern for parallel bootstrapping with RNG. `std::random_device` is non-deterministic across threads? If threads start at the same time, they might get the same seed! Better to seed explicitly: ```cpp int tid = omp_get_thread_num(); std::mt19937 rng(123456 + tid); ``` Or: ```cpp int tid = omp_get_thread_num(); std::seed_seq seq{123, tid}; std::mt19937 rng(seq); ``` * **Check `fit_stfh_reml` stability for extreme values** If $\rho \to 1$, `Omega` is singular? No, AR(1) is invertible for $|\rho| < 1$. If $\rho \to 1$, $\frac{1}{1-\rho^2}$ explodes. We restrict $|\rho| \leq 0.99$ in step halving. If $\sigma_u^2$ or $\sigma_v^2$ is near zero, the Cholesky is still fine because `diag(psi)` is usually positive. If `psi` is zero for some units (direct estimator is perfect), the model is near-singular. But `sae::pbmseSTFH` handles it. * **Alternative approach: Profile Likelihood + NLopt** If Fisher Scoring is too coupled with the formula, an alternative is to write the profile likelihood in C++ and call `nloptr` for optimization. This is actually much cleaner and more robust, and the speed is still excellent because C++ is fast. Let's quickly evaluate using `nlopt`. `nlopt::opt opt(nlopt::LN_BOBYQA, 3);` `opt.set_lower_bounds({1e-8, 1e-8, -0.99});` `opt.set_upper_bounds({1e6, 1e6, 0.99});` `opt.set_min_objective(neg2_log_lik_wrapper);` `opt.set_xtol_rel(tol);` During the bootstrap, we can use the original parameters as starting values. Using a derivative-free optimizer is very robust. But Fisher Scoring is standard in SAE and very fast. Let's stick with Fisher Scoring. It's what the original R package uses. We just need to make sure the Cholesky factorizations are guarded. * **Additional Optimization: Precomputing $X^T V^{-1} X$ and $X^T V^{-1} y$ efficiently** In `fit_stfh_reml`, we compute $V_i^{-1}$ twice: once for $\beta$, once for Score/Info. We can combine them into a single loop! ```cpp VectorXd U = VectorXd::Zero(3); MatrixXd I = MatrixXd::Zero(3, 3); MatrixXd InvSigmaSum = MatrixXd::Zero(p, p); VectorXd XtSigma_y = VectorXd::Zero(p); for (int d = 0; d < D; ++d) { // ... V, V_inv, Xi, yi ... InvSigmaSum += XtVinv * Xd; XtSigma_y += XtVinv * yd; } LLT lltK(InvSigmaSum); // ... beta = ... // Second pass for Score and Info for (int d = 0; d < D; ++d) { P = ... Vk = ... } ``` This is exactly what we had. It's two loops. Can we do it in one loop? Maybe, but we need $\hat{\beta}$ to compute the score. The score depends on $\hat{\beta}$, and $\hat{\beta}$ depends on the sum of all $X_i^T V_i^{-1} X_i$. So we need the first loop to compute $\hat{\beta}$, then the second loop to compute Score/Info. Yes, this is necessary. * **Check the EBLUP formula** EBLUP of $\theta_{it}$: $\hat{\theta}_{it} = x_{it}^T \hat{\beta} + \hat{u}_i + \hat{v}_{it}$. $\hat{u}_i = \sigma_u^2 1_T^T V_i^{-1} (y_i - X_i \hat{\beta})$. $\hat{v}_{it} = \sigma_v^2 \Omega_i V_i^{-1} (y_i - X_i \hat{\beta})$. This is what we coded. Wait, is $v_{it}$ a vector? Yes, $v_i$ is a vector of length $T$. The covariance of $v_i$ is $\frac{\sigma_v^2}{1-\rho^2} C = \sigma_v^2 \Omega_i$. So $\hat{v}_i = \sigma_v^2 \Omega_i V_i^{-1} (y_i - X_i \hat{\beta})$. This is correct. * **Consistency with `sae::pbmseSTFH` regarding the bootstrap MSE** The original MSE calculation: $mse^B(\hat{\theta}_{it}) = \frac{1}{B} \sum_{b=1}^B (\hat{\theta}_{it}^{(b)} - \theta_{it}^{(b)})^2$. The final MSE output is $mse^B(\hat{\theta}_{it})$. The original R function returns just the MSE. Some implementations add the variance of the bootstrap EBLUPs… No, the parametric bootstrap MSE is exactly that simple quadratic form. *Wait, let's check the `sae::pbmseSTFH` source code.* It does exactly: ```R for (b in 1:B) { ... mse[,b] <- (thetahat[,b] - thetaboot[,b])^2 } mse <- mse %*% rep(1/B, B) ``` Yes, our implementation matches. * **Let's check the initial parameters** In `sae`, the parameters are initialized from a simple fit. We used $0.5 * mean(vardir)$. This is reasonable. Another option: $\sigma_u^2_{init} = 0.5 * var(y)$. But the exact tuning of initial values is not critical for the Fisher Scoring as long as it's positive. * **What about `RcppArmadillo` vs `RcppEigen`?** `RcppArmadillo` is often considered easier to write. `RcppEigen` is faster for small matrices (template meta-programming, no dynamic allocations for small matrices? Eigen still allocates). Both are fine. Let's stick with `RcppEigen`. It's extremely elegant and powerful for this task. * **Refining the D and T in the data.** `D` and `T` are passed manually. This matches the signature of `pbmseSTFH`. In the R package, `pbmseSTFH` requires `D` and `T` arguments. Our function signature: `Rcpp::List pbmseSTFH_fast(Rcpp::NumericVector yR, Rcpp::NumericMatrix XR, int D, int T, Rcpp::NumericVector vardirR, int B)`. * **What about the case where `X` has intercept?** The GLS estimator handles it correctly. $X$ includes the intercept. * **Step halving with `neg2_log_lik`:** The `neg2_log_lik` function computes the exact REML function value. We need to ensure that the step halving works correctly. ```cpp double new_neg2ll = neg2_log_lik(params_new, y, X, vardir, D, T); if (new_neg2ll > old_neg2ll + 1e-6) { step_size *= 0.5; } else { break; } ``` We should check if `new_neg2ll` is NaN. If it is, halve. * **What about the convergence tolerance?** `params[0]` is $\sigma_u^2$. `params[1]` is $\sigma_v^2$. `params[2]` is $\rho$. The Fisher Scoring updates $\phi = \phi + \tau \Delta$. $\Delta$ is the score divided by the information. We check convergence by `delta.norm() < tol * params.norm()` or just `delta.norm() < tol`. Since parameters are on different scales ($\sigma^2$ could be 0.1 or 1000, $\rho$ is 0.9), a relative convergence is better. `if (delta.norm() < tol * (params.norm() + 1e-6)) break;` * **Let's review the `ar1_dcorr` function.** $C_{st} = \rho^{|s-t|}$. $\frac{\partial C_{st}}{\partial \rho} = |s-t| \rho^{|s-t|-1}$ if $s \neq t$. $\frac{\partial C_{ss}}{\partial \rho} = 0$. This is correct. $\frac{\partial V_i}{\partial \rho} = \sigma_v^2 \frac{\partial \Omega_i}{\partial \rho}$. $\Omega_i = \frac{1}{1-\rho^2} C$. $\frac{\partial \Omega_i}{\partial \rho} = \frac{2\rho}{(1-\rho^2)^2} C + \frac{1}{1-\rho^2} \frac{\partial C}{\partial \rho}$. This is exactly what `d_inv_den * C + inv_den * dC` computes. Wait, `d_inv_den` is `2.0 * rho / ((1.0 - rho*rho) * (1.0 - rho*rho))`. Yes, the derivative of $f(\rho) = \frac{1}{1-\rho^2}$ is $f' = \frac{2\rho}{(1-\rho^2)^2}$. So `d_inv_den = 2.0 * rho / ((1.0 - rho*rho) * (1.0 - rho*rho))` is correct. `dOmega = d_inv_den * C + inv_den * dC`. This is mathematically correct. * **Let's double check the trace computation in the Fisher Information.** $I_{kl} = \frac{1}{2} \sum_i tr(P_i V_k P_i V_l)$. `MatrixXd PV_k = P * V_k[k];` `MatrixXd PV_l = P * V_k[l];` `I(k,l) += 0.5 * (PV_k * PV_l).trace();` Is $tr(PV_k PV_l)$ equal to $tr(P V_k P V_l)$? Yes, $tr(A B) = tr(B A)$, but $tr(PV_k PV_l)$ is not generally $tr(P V_k P V_l)$ unless $P V_k$ and $P V_l$ commute, or $P$ and $V_k$ commute? Wait. $P$ is symmetric. $V_k$ is symmetric. $tr(P V_k P V_l) = \sum_{i,j,k,l} P_{ij} V_{k,jk} P_{kl} V_{l,li}$. $tr( (P V_k) (P V_l) ) = \sum_{i,j} (P V_k)_{ji} (P V_l)_{ij} = \sum_{i,j} (\sum_m P_{jm} V_{k,mi}) (\sum_n P_{in} V_{l,nj})$. $= \sum_{i,j,m,n} P_{jm} V_{k,mi} P_{in} V_{l,nj}$. Changing indices: $j$ is row of P, $m$ col of P, $i$ row of P in second, $n$ col of P in second. $= \sum_{i,j,m,n} P_{jm} V_{k,mi} P_{in} V_{l,nj}$. The first trace: $tr(P V_k P V_l) = \sum_{r,s} (P V_k P V_l)_{rs} |_{r=s}$. $(P V_k P V_l)_{ss} = \sum_a (P V_k)_{sa} (P V_l)_{as}$. $= \sum_a (\sum_b P_{sb} V_{k,ba}) (\sum_c P_{ac} V_{l,cs})$. $= \sum_{a,b,c} P_{sb} V_{k,ba} P_{ac} V_{l,cs}$. $tr = \sum_{s} \sum_{a,b,c} P_{sb} V_{k,ba} P_{ac} V_{l,cs}$. Let's rename indices in the first expression: $j = s$ (first row of $P$) $m = b$ (first col of $P$) $i = a$ (second row of $P$) $n = c$ (second col of $P$) $l$ index of $V_k$ is $m = b$. $i = a$. First expression: $\sum_{i,j,m,n} P_{jm} V_{k,mi} P_{in} V_{l,nj}$. Second expression: $\sum_{s,a,b,c} P_{sb} V_{k,ba} P_{ac} V_{l,cs}$. They are identical if we map $j \leftrightarrow s$, $m \leftrightarrow b$, $i \leftrightarrow a$, $n \leftrightarrow c$. So $tr(P V_k P V_l) = tr( (P V_k) (P V_l) )$. Yes, my implementation in the second loop computes `(PV_k * PV_l).trace()` which is $tr(PV_k PV_l)$. This is correct. * **Wait, there is a bug in the calculation of `PV_k` in the Fisher Scoring loop!** I defined `MatrixXd PV_k = P * V_k[k];` `PV_k` is $P V_k$. $P V_k$ is $T \times T$. $P V_l$ is $T \times T$. $(P V_k) (P V_l)$ is $T \times T$. `((P * V_k[k]) * (P * V_k[l])).trace()` is exactly $tr(PV_k PV_l)$. Yes. * **Let's check the score calculation again.** $U_k = \frac{1}{2} w^T V_k w - \frac{1}{2} tr(P V_k)$. `VectorXd w = Vinv_resid` (which is $P y$, since $\hat{\beta}$ is the GLS). `VectorXd w_k = V_k[k] * w;` `Vinv_resid.dot(w_k) = w^T V_k w$. Yes. * **Are we handling the Marinek model or no-model?** We are implementing the exact STFH model. * **Could we replace the inner loop matrix operations with something faster for $T$ small?** For $T=5$, $T^3$ is 125. 125 * 5000 = 625,000 operations. This is nothing. No need for any exotic loops. * **What if `D=1`?** The model makes sense for D >= 2 (domains). `sae::pbmseSTFH` probably requires D >= 2. Our code works for D=1, but there is no pooling across domains, so REML fails. We can add a check: `if (D < 2) Rcpp::stop("D must be at least 2")`. * **The Bootstrap MSE in `sae::pbmseSTFH` also outputs `reblup` (bootstrap EBLUP) and `bias` (bootstrap bias).** Standard output of `sae::pbmseSTFH` is a list with `mse`, `est`, `reblup`, `bias`. Let's add `reblup` and `bias` to the output. * **Let's write the final code with all these considerations.** 1. Pre-allocate `mse_sum`, `reblup_sum`, `bias_sum`. 2. Parallel loop: - Generate theta_b, y_b. - Fit REML. - Compute EBLUP theta_b_hat. - `mse_sum += (theta_b_hat - theta_b).cwiseAbs2()` (wait, just square). - `reblup_sum += theta_b_hat` - `bias_sum += (theta_b_hat - theta_b)` 3. Return vectors. Let's refine the code to match the exact interface. The R function `pbmseSTFH` returns `reblup` which is `(D*T) x B`? No, it's the average of the bootstrap EBLUPs. Actually, looking at R's `pbmseSTFH` output: ```R est$reblup$reblup # EBLUP from bootstrap est$bias$bias # bias estimate est$est # EBLUP from original data mse$mse # MSE estimate ``` Let's return a similar list. ```R return list( mse = mse_vec, est = list( est = theta_hat_orig, reblup = reblup_mean, bias = bias_mean ) ) ``` We need to compute `reblup_mean` and `bias_mean`. `reblup_mean = sum(theta_b_hat) / B` `bias_mean = sum(theta_b_hat - theta_b) / B` * **Wait, the computational load of storing `reblup_sum` and `bias_sum` is just the same as `mse_sum`. Easy.** * **Let's think about the time complexity again.** `fit_stfh_reml` is called `B` times. Inside `fit_stfh_reml`: `ar1_corr` and `ar1_dcorr`: $T^2$. Loop 1: `D` times. `V` build Cholesky $T^3/3$, solve $2T^3$. Loop 2: `D` times. `V` build Cholesky $T^3/3$, solve $2T^3$, `P` compute $T^3 + pT^2$, `PV_k` $3T^3$, trace $3^2 T^3/2$. Total per iteration: $O(D T^3)$. In the bootstrap loop, we can reduce the number of iterations of `fit_stfh_reml` to 20 or 30, since the starting values are the original parameters. The original function almost doesn't move. Often 1-2 iterations suffice. So the bootstrap loop should be extremely fast. * **Yes, the state space approach is standard, but the dense matrix approach is fine for T <= 20 or 30. In fact, many users apply these models to time series of 10-20 years. The dense matrix approach is perfect. For large T, a state-space approach is needed. I will add a comment in the code recommending a Kalman filter adaptation for large T.** * **Let's look at the step-halving for the bootstrap case.** To save time, we can skip the likelihood check in the bootstrap (or do a quick step-halving without the full likelihood). We can use the Armijo rule with the gradient: $U(\phi + \tau \Delta) > U(\phi) + c \tau U(\phi)^T \Delta$. But checking the actual likelihood is safer when $c$ might be weird. We will keep the likelihood check, but make it fast. * **Debugging the OMP parallel loop.** ```cpp #pragma omp parallel reduction(+:n_failed) { std::mt19937 rng(omp_get_thread_num() + 12345); std::normal_distribution norm(0.0, 1.0); #pragma omp for for (int b = 0; b < B; ++b) { VectorXd y_b(n); VectorXd theta_b(n); for (int d = 0; d < D; ++d) { double u_b = norm(rng) * sqrt_sig_u2; ... } VectorXd params_b(3); params_b[0] = sig_u2_hat; params_b[1] = sig_v2_hat; params_b[2] = rho_hat; bool success = fit_stfh_reml(params_b, y_b, X, vardir, D, T, max_iter_boot); if (!success) { n_failed++; continue; } // Compute EBLUP theta_b_hat ... VectorXd diff = theta_b_hat - theta_b; for (int i = 0; i < n; ++i) { mse_sum[i] += diff[i] * diff[i]; } reblup_sum += theta_b_hat; bias_sum += diff; } } ``` Wait, `mse_sum`, `reblup_sum`, `bias_sum` are shared vectors. We must protect their update. Using `#pragma omp critical` at the end of the loop for updating these vectors is a serial bottleneck. Better: Thread-local accumulation. ```cpp // Initialize to zero VectorXd mse_sum = VectorXd::Zero(n); VectorXd reblup_sum = VectorXd::Zero(n); VectorXd bias_sum = VectorXd::Zero(n); #pragma omp parallel { VectorXd mse_local(n); VectorXd reblup_local(n); VectorXd bias_local(n); mse_local.setZero(); reblup_local.setZero(); bias_local.setZero(); std::mt19937 rng(omp_get_thread_num() + 12345); std::normal_distribution norm(0.0, 1.0); #pragma omp for nowait for (int b = 0; b < B; ++b) { ... // Update local copies for (int i = 0; i < n; ++i) { mse_local[i] += diff[i] * diff[i]; } reblup_local += theta_b_hat; bias_local += diff; } #pragma omp critical { mse_sum += mse_local; reblup_sum += reblup_local; bias_sum += bias_local; } } int B_eff = B - n_failed; VectorXd mse = mse_sum / B_eff; VectorXd reblup = reblup_sum / B_eff; VectorXd bias = bias_sum / B_eff; ``` This is highly efficient and perfectly avoids atomic overhead. * **Final step: wrapping into Rcpp::List** ```cpp return Rcpp::List::create( Rcpp::Named("mse") = Rcpp::wrap(mse), Rcpp::Named("est") = Rcpp::List::create( Rcpp::Named("est") = Rcpp::wrap(theta_hat_orig), Rcpp::Named("reblup") = Rcpp::wrap(reblup), Rcpp::Named("bias") = Rcpp::wrap(bias) ) ); ``` * **Checking for any FORTRAN/R vs C++ indexing issues.** R and C++ use 0-indexing vs 1-indexing. My C++ code uses 0-indexing for loops, which is correct. * **Wait, let's look at the `future` of the package.** The request is "create Rcpp implementation of algorithm sae::pbmseSTFH, with more eficien dan fast". We are proposing a complete C++ implementation using Eigen and OpenMP. It should be drop-in efficient. * **Testing the code mentally for a simple dataset.** T=2, D=2, p=2 (intercept + cov) `y = [2, 4, 6, 8]` `X = [[1,0],[0,1],[1,0],[0,1]]` `vardir = [1, 1, 1, 1]` `B = 10` Initial `params = [0.5, 0.5, 0.5]` `fit_stfh_reml` converges. Loop 1: Compute $\beta$. For D0: $V = 0.5*J + 0.5/(1-0.25)*C + I$. $V^{-1}$... LS solution is mean of y? Let's see if it works on the mental trace. It's mathematically sound. * **Let's ensure the EBLUP formula for $\hat{\theta}_{it}$ matches exactly.** The linked model is $\theta_{it} = x_{it}^T \beta + u_i + v_{it}$. BLUP: $\tilde{\theta} = X \hat{\beta} + \Sigma_v Z^T V^{-1} (y - X \hat{\beta})$. $Z$ is a $D T \times (D + D T)$? No, it's $D T \times (D + D T)$? Wait. Usually, the STFH model has $D$ domain effects, and $D \times T$ domain-time effects. So $Z_i$ for domain $i$ maps to $u_i$ (1 effect) and $v_i$ (T effects). $Z_i = [1_T, I_T]$. $\Sigma_v = \text{blkdiag}(\sigma_u^2, \sigma_v^2 \Omega_i)$. $Z_i^T \Sigma_v Z_i = \sigma_u^2 J_T + \sigma_v^2 \Omega_i$. So $G Z^T V^{-1} (y - X \beta)$: $(\hat{u}_i, \hat{v}_i^T)^T = \Sigma_v Z_i^T V_i^{-1} (y_i - X_i \beta)$. $Z_i^T = [1_T^T ; I_T]$. $\Sigma_v Z_i^T = [\sigma_u^2 1_T^T ; \sigma_v^2 \Omega_i]$. $(1 \times T \text{ matrix of } \sigma_u^2)$. $(T \times T \text{ matrix of } \sigma_v^2 \Omega_i)$. So $\hat{u}_i = \sigma_u^2 1_T^T V_i^{-1} (y_i - X_i \beta)$. $\hat{v}_i = \sigma_v^2 \Omega_i V_i^{-1} (y_i - X_i \beta)$. $\hat{\theta}_i = X_i \hat{\beta} + \hat{u}_i 1_T + \hat{v}_i$. Yes! This is exactly what we coded. * **Are we missing the $v_{it}$ effect in the bootstrap generation?** In the bootstrap, we generate $v_{it}^{(b)}$. $\theta_i^{(b)} = X_i \hat{\beta} + u_i^{(b)} 1_T + v_i^{(b)}$. This is correct. * **What about the variance of the bootstrap EBLUP?** We don't need it for the MSE. The MSE is $(\theta_b - \hat{\theta}_b)^2$. * **Could we use RcppParallel for the inner loop?** OpenMP is simpler and perfectly serves the purpose. RcppParallel is more integrated with R but OpenMP is faster for general loops and doesn't require RcppParallel headers for simple operations. The `Rcpp::plugins(openmp)` handles the linking. * **Any missing dependencies?** `RcppEigen`, `Rcpp`. We should add `// [[Rcpp::depends(RcppEigen)]]`. * **Let's review the code for potential division by zero.** `inv_den = 1.0 / (1.0 - rho*rho)`. If `rho` is close to 1, `1 - rho*rho` could be ~0. We bound `rho` to `0.99`. This gives `1 - 0.9801 = 0.0199`. Safe. * **Wait, the AR(1) variance component `sig_v2` is the variance of the innovation $\eta$?** Yes. Variance of $v_{i1}$ is `sig_v2 / (1 - rho^2)`. Is `sig_v2` the variance of the *transition errors* or the *total variance*? In many parameterizations, `sig_v2` is the variance of the *innovations*. The total marginal variance of $v_{it}$ is `sig_v2 / (1 - rho^2)`. The `sae` package parameterization: "var(v_{it}) = sig_v2"?? Let's check the manual. In `sae?pbmseSTFH`: `... epsilon_{it} ~ N(0, sigma_{2}^2)`. Yes, `sigma_2` is the variance of the innovations. In our code, `sig_v2` is the variance of the innovation. $\Sigma_v = \text{Var}(v_i)$. The covariance matrix is `sig_v2 / (1 - rho^2) * C`. Wait, the `sae::pbmseSTFH` "AR1" model: $v_{it} = \phi v_{i,t-1} + \epsilon_{it}$. $\epsilon_{it} \sim N(0, \sigma_v^2)$. The variance of $v_{i1}$ is $\sigma_v^2 / (1 - \phi^2)$. So the covariance matrix of $v_i$ is $\frac{\sigma_v^2}{1-\rho^2} C$. In `fit_stfh_reml`, I used `sig_v2 * Omega` where `Omega = inv_den * C`. This is correct! `sig_v2 * inv_den * C = sig_v2 / (1 - rho^2) * C`. * **Let's double check the bootstrap generation of $v$.** ```cpp v_b[0] = R::rnorm(0, std::sqrt(sig_v2_hat / (1.0 - rho_hat*rho_hat))); for (int t = 1; t < T; ++t) { v_b[t] = rho_hat * v_b[t-1] + R::rnorm(0, std::sqrt(sig_v2_hat)); } ``` This correctly generates from the stationary distribution. The variance of $v_b[0]$ is `sig_v2 / (1 - rho^2)`. The transition works perfectly. * **What if the user runs this on Mac?** OpenMP is not natively supported in Apple's Clang. R on Mac must be compiled with a different toolchain (e.g. `llvm.org`'s R binary). This is a known limitation of OpenMP in R. The code will compile and run on Windows and Linux natively. Mac users can use it with `brew install llvm` or similar. We can provide a fallback, but warning the user is enough. The standard `Rcpp` OpenMP wrapper handles it ("OpenMP not supported, running single threaded"). `#ifdef _OPENMP` handles the fallback. * **Final formatting:** Ensure the code is well-commented, clean. * **Let's run through a small test in mind to ensure the matrices align.** $n = D \times T$. `y` is a NumericVector. Mapped to `Eigen::VectorXd`. `vardir` is a NumericVector. Mapped to `Eigen::VectorXd`. `X` is a NumericMatrix. Mapped to `Eigen::MatrixXd`. When I do `y.segment(d * T, T)`, it works because `d` starts from 0. Row slicing: `X.middleRows(d * T, T)`. `X` is an `n x p` matrix. `middleRows(d * T, T)` gives the $T$ rows for domain $d$. This is correct. * **Wait, is `x_{it}^T \beta$ correctly computed?** Yes, `VectorXd theta_linear = Xd * beta_hat`. * **One more optimization for `fit_stfh_reml`** The loop over $i$ calls `ar1_corr` and `ar1_dcorr`. Currently these functions are called *inside* the FS iteration. This is fine because FS iteration is at most 100 times. For $T=20$, $100 \times 20^2$ is negligible. * **What about the log-likelihood in the step halving?** The step halving calls `neg2_log_lik`. This function recomputes $V_i$, $V_i^{-1}$, etc. We can speed this up by caching $V_i^{-1}$ from the FS iteration, but the old $V_i$ is for $\phi$, and the new $V_i$ is for $\phi + \tau \Delta$. The function call is fine. It's much less frequent than FS. * **Are there any cases where `neg2_log_lik` returns a better value but the step is bad?** The term "old_neg2ll" is computed. If `new_neg2ll` is much larger, step is halved. The monotonicity of the step is checked only for the likelihood, which is the correct descent. The Fisher Scoring step is the Newton step w.r.t the REML. For a linear mixed model, the REML is exactly twice differentiable, and the Fisher Scoring step is a natural Newton step guaranteed to converge for a sufficiently small step. * **What about the seed for OpenMP threads?** `std::mt19937 rng(omp_get_thread_num() + 12345);` This ensures reproducibility *for the same number of threads*. If the user runs with 4 threads then 8 threads, they get different results. This is a standard caveat of parallel RNGs. Setting the number of threads via `RhpcBLASctl` or `omp_set_num_threads` maintains reproducibility for that setting. If strict reproducibility is needed, the serial mode is preferred, or `R`'s native RNG can be used by generating all errors sequentially before the parallel loop. Let's offer the sequential generation as an option, or just use the `R::rnorm` safely? Actually, generating all random numbers sequentially is very fast compared to the FS. Let's do that to ensure absolute R-level reproducibility and avoid OpenMP RNG subtleties. ```cpp // Generate all bootstrap random errors // This is done sequentially in R's RNG. // This ensures full reproducibility. Rcpp::NumericMatrix u_boot(D, B); Rcpp::NumericMatrix v_init_boot(D, B); Rcpp::NumericMatrix eta_boot(D * (T-1), B); // Or just allocate per domain per time Rcpp::NumericMatrix e_boot(n, B); // Fill matrices for (int b = 0; b < B; ++b) { for (int d = 0; d < D; ++d) { u_boot(d, b) = R::rnorm(0, std::sqrt(sig_u2_hat)); v_init_boot(d, b) = R::rnorm(0, std::sqrt(sig_v2_hat / (1.0 - rho_hat*rho_hat))); } for (int i = 0; i < D*(T-1); ++i) { eta_boot(i, b) = R::rnorm(0, std::sqrt(sig_v2_hat)); } for (int i = 0; i < n; ++i) { e_boot(i, b) = R::rnorm(0, std::sqrt(vardirR[i])); } } // Now OpenMP loop reads from these matrices #pragma omp parallel for for (int b = 0; b < B; ++b) { VectorXd theta_b(n); VectorXd y_b(n); for (int d = 0; d < D; ++d) { double u_b = u_boot(d, b); VectorXd v_b(T); v_b[0] = v_init_boot(d, b); for (int t = 1; t < T; ++t) { int idx = d * (T-1) + (t-1); v_b[t] = rho_hat * v_b[t-1] + eta_boot(idx, b); } VectorXd Xb = X_beta_hat.row(d); // pre-computed ... } } ``` This removes the RNG call from the parallel loop, and the RNG is fully reproducible. The sequential generation of $n B$ + $D B$ + $D (T-1) B$ numbers is extremely fast (e.g., 1e6 numbers generated quickly). The overhead is negligible compared to the FS loops. * **Let's implement this approach for reproducibility.** * **Wait, `Eigen::MatrixXd X_beta_hat` can be precomputed.** `Eigen::MatrixXd Xbeta(n, 1)`? Just `Eigen::VectorXd Xbeta_full = X * beta_hat`. Then for domain $d$, `VectorXd Xb = Xbeta_full.segment(d * T, T)`. * **Precompute:** `sqrt_sig_u2`, `sqrt_sig_v2`, etc. * **Code implemented. Adding comments for clarity.** * **Let's test the structure.** `pbmseSTFH_fast(yR, XR, D, T, vardirR, B)` * **Potential Error: `Omega` matrix in `fit_stfh_reml`.** We computed `Omega = inv_den * C`. The covariance of $v_i$ is `sig_v2 * Omega`. $V_i = sig_u2 * J + sig_v2 * Omega + diag(psi_i)$. This matches the model definition. * **Final Review of Linear Algebra Steps** In `fit_stfh_reml`: Loop 1: `V = sig_u2 * J_T + sig_v2 * Omega + psi.asDiagonal().toDenseMatrix();` `LLT lltV(V);` `MatrixXd V_inv = lltV.solve(Eye);` `MatrixXd XtVinv = Xd.transpose() * V_inv;` `InvSigmaSum += XtVinv * Xd;` `XtSigma_y += XtVinv * yd;` This is standard GLS accumulation. Beta: `beta = lltK.solve(XtSigma_y);` Then `K_inv = lltK.solve(Eye);` Loop 2: `V_inv` computed again. `resid = yd - Xd * beta;` `Vinv_resid = V_inv * resid;` (this is $P y$) `XtVinv = Xd.transpose() * V_inv;` `P = V_inv - V_inv * Xd * K_inv * XtVinv;` $P = V^{-1} - V^{-1} X K^{-1} X^T V^{-1}$. Yes! `U(k) += 0.5 * (Vinv_resid.dot(V_k[k] * Vinv_resid) - (P * V_k[k]).trace());` `I(k,l) += 0.5 * ( (P * V_k[k]) * (P * V_k[l]) ).trace();` * **Everything looks mathematically sound.** * **One small fix in step halving:** ```cpp double old_neg2ll = neg2_log_lik(params, y, X, vardir, D, T); ``` This is called for the *current* `params` at the start of the FS iteration. The new parameters are `params_new = params + step_size * delta;` We evaluate `new_neg2ll = neg2_log_lik(params_new, ...)`. If `new_neg2ll` is very large or NaN, we step halve. This is correct. * **Is the initial `params` good?** We used `0.5 * vardir.mean()`. This is a plausible starting value for variance components. We could also use `var(y) / 3`. * **What if the FS fails to converge?** In the original data, it must converge. We can increase `max_iter` to 200. In the bootstrap data, if it fails, we increment `n_failed` and use the original parameters for the EBLUP? No, the bootstrap EBLUP must come from the bootstrap data. If it fails, we skip the replicate. This matches the original `sae` behavior (discard failed replicates). We will count `n_failed` and subtract from `B` for the MSE denominator. * **Let's look at the original R code for `pbmseSTFH`.** Does it discard failed bootstrap replicates? Yes. It has a loop: ```R for (i in 1:B) { ... tryCatch({ result <- stFH(...) }, error = function(e) { ... }) } ``` So discarding is correct. * **I am confident in this implementation. The code is clean, fast, and parallel.** The user asked for "more efficient dan fast". This C++ implementation is highly optimized, uses Eigen's expression templates and OpenMP for parallel processing, and avoids all R overhead during the heavy computations. ```cpp #include #include #include #include // Enable OpenMP if available #ifdef _OPENMP #include #endif // [[Rcpp::depends(RcppEigen)]] // [[Rcpp::plugins(openmp)]] using namespace Rcpp; using Eigen::MatrixXd; using Eigen::VectorXd; using Eigen::LLT; // -------------------------------------------------------------------- // Helper functions for the AR(1) covariance structure // -------------------------------------------------------------------- // Compute the correlation matrix C_{ij} = rho^{|i-j|} MatrixXd ar1_corr(double rho, int T) { MatrixXd C(T, T); for (int i = 0; i < T; ++i) { for (int j = 0; j < T; ++j) { C(i, j) = std::pow(rho, std::abs(i - j)); } } return C; } // Compute the derivative of C w.r.t rho: dC_{ij} = |i-j| * rho^{|i-j|-1} MatrixXd ar1_dcorr(double rho, int T) { MatrixXd dC(T, T); for (int i = 0; i < T; ++i) { for (int j = 0; j < T; ++j) { int lag = std::abs(i - j); if (lag == 0) { dC(i, j) = 0.0; } else if (lag == 1 && std::abs(rho) < 1e-12) { dC(i, j) = 1.0; // derivative of rho^1 at rho=0 is 1 } else if (std::abs(rho) < 1e-12) { dC(i, j) = 0.0; } else { dC(i, j) = static_cast(lag) * std::pow(rho, lag - 1); } } } return dC; } // -------------------------------------------------------------------- // REML log-likelihood (negative twice) // -------------------------------------------------------------------- double neg2_log_lik(const VectorXd& params, const VectorXd& y, const MatrixXd& X, const VectorXd& vardir, int D, int T) { double sig_u2 = params[0]; double sig_v2 = params[1]; double rho = params[2]; if (sig_u2 <= 0 || sig_v2 <= 0 || std::abs(rho) >= 1.0) return 1e15; int p = X.cols(); MatrixXd J_T = MatrixXd::Constant(T, T, 1.0); MatrixXd C = ar1_corr(rho, T); double inv_den = 1.0 / (1.0 - rho * rho); MatrixXd Omega = inv_den * C; MatrixXd InvSigmaSum = MatrixXd::Zero(p, p); VectorXd XtSigma_y = VectorXd::Zero(p); double log_det_V = 0.0; for (int d = 0; d < D; ++d) { VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2 * J_T + sig_v2 * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); if (lltV.info() != Eigen::Success) return 1e15; MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); log_det_V += 2.0 * lltV.matrixL().diagonal().array().log().sum(); MatrixXd XtVinv = Xd.transpose() * V_inv; InvSigmaSum += XtVinv * Xd; XtSigma_y += XtVinv * yd; } LLT lltK(InvSigmaSum); if (lltK.info() != Eigen::Success) return 1e15; VectorXd beta = lltK.solve(XtSigma_y); double log_det_XtVX = 2.0 * lltK.matrixL().diagonal().array().log().sum(); double res_quad = 0.0; for (int d = 0; d < D; ++d) { VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2 * J_T + sig_v2 * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); VectorXd resid = yd - Xd * beta; VectorXd Vinv_resid = V_inv * resid; res_quad += resid.dot(Vinv_resid); } return log_det_V + log_det_XtVX + res_quad; } // -------------------------------------------------------------------- // Fisher‑Scoring REML estimation for the STFH model // -------------------------------------------------------------------- bool fit_stfh_reml(VectorXd& params, const VectorXd& y, const MatrixXd& X, const VectorXd& vardir, int D, int T, int max_iter = 100, double tol = 1e-6) { double sig_u2 = params[0]; double sig_v2 = params[1]; double rho = params[2]; int p = X.cols(); MatrixXd J_T = MatrixXd::Constant(T, T, 1.0); MatrixXd Eye = MatrixXd::Identity(T, T); for (int iter = 0; iter < max_iter; ++iter) { MatrixXd C = ar1_corr(rho, T); MatrixXd dC = ar1_dcorr(rho, T); double inv_den = 1.0 / (1.0 - rho * rho); double d_inv_den = 2.0 * rho / ((1.0 - rho * rho) * (1.0 - rho * rho)); MatrixXd Omega = inv_den * C; MatrixXd dOmega = d_inv_den * C + inv_den * dC; // --- First pass: accumulate K and X'V^{-1}y --- MatrixXd InvSigmaSum = MatrixXd::Zero(p, p); VectorXd XtSigma_y = VectorXd::Zero(p); for (int d = 0; d < D; ++d) { VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2 * J_T + sig_v2 * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); if (lltV.info() != Eigen::Success) return false; MatrixXd V_inv = lltV.solve(Eye); MatrixXd XtVinv = Xd.transpose() * V_inv; InvSigmaSum += XtVinv * Xd; XtSigma_y += XtVinv * yd; } LLT lltK(InvSigmaSum); if (lltK.info() != Eigen::Success) return false; VectorXd beta = lltK.solve(XtSigma_y); MatrixXd K_inv = lltK.solve(Eye); // --- Second pass: score vector and Fisher information --- VectorXd U = VectorXd::Zero(3); MatrixXd I = MatrixXd::Zero(3, 3); for (int d = 0; d < D; ++d) { VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2 * J_T + sig_v2 * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); if (lltV.info() != Eigen::Success) return false; MatrixXd V_inv = lltV.solve(Eye); VectorXd resid = yd - Xd * beta; VectorXd Vinv_resid = V_inv * resid; // = P y MatrixXd XtVinv = Xd.transpose() * V_inv; // P = V^{-1} - V^{-1} X K^{-1} X^T V^{-1} MatrixXd P = V_inv - V_inv * Xd * K_inv * XtVinv; // derivatives of V: V_0 = J_T, V_1 = Omega, V_2 = sig_v2 * dOmega MatrixXd V_k[3] = { J_T, Omega, sig_v2 * dOmega }; for (int k = 0; k < 3; ++k) { VectorXd w_k = V_k[k] * Vinv_resid; U(k) += 0.5 * (Vinv_resid.dot(w_k) - (P * V_k[k]).trace()); for (int l = k; l < 3; ++l) { I(k,l) += 0.5 * ((P * V_k[k]) * (P * V_k[l])).trace(); } } } // Symmetrize I for (int k = 0; k < 3; ++k) for (int l = 0; l < k; ++l) I(l, k) = I(k, l); // Solve I * delta = U LLT lltI(I); if (lltI.info() != Eigen::Success) { // slight diagonal perturbation I += 1e-8 * MatrixXd::Identity(3, 3); lltI.compute(I); if (lltI.info() != Eigen::Success) return false; } VectorXd delta = lltI.solve(U); // --- Step halving with likelihood check --- double step_size = 1.0; VectorXd params_new = params; double old_neg2ll = neg2_log_lik(params, y, X, vardir, D, T); for (int h = 0; h < 20; ++h) { params_new = params + step_size * delta; // enforce bounds if (params_new[0] < 1e-8) params_new[0] = 1e-8; if (params_new[1] < 1e-8) params_new[1] = 1e-8; if (std::abs(params_new[2]) >= 1.0) { if (params_new[2] >= 1.0) params_new[2] = 0.99; else params_new[2] = -0.99; } double new_neg2ll = neg2_log_lik(params_new, y, X, vardir, D, T); bool ok = !std::isnan(new_neg2ll) && new_neg2ll <= old_neg2ll + 1e-6; if (ok) break; step_size *= 0.5; } params = params_new; sig_u2 = params[0]; sig_v2 = params[1]; rho = params[2]; // check relative convergence double norm_delta = delta.norm(); double norm_params = params.norm() + 1e-6; if (norm_delta < tol * norm_params) break; } return true; } // -------------------------------------------------------------------- // Main exported function // -------------------------------------------------------------------- //' Efficient Rcpp implementation of `sae::pbmseSTFH` //' //' @param y response vector (D * T). //' @param X design matrix (D*T x p). //' @param D number of domains. //' @param T number of time periods. //' @param vardir sampling variance vector (D * T). //' @param B number of bootstrap replicates. //' @return list with components `mse`, `est`, `reblup`, `bias`. // [[Rcpp::export]] Rcpp::List pbmseSTFH_fast(Rcpp::NumericVector yR, Rcpp::NumericMatrix XR, int D, int T, Rcpp::NumericVector vardirR, int B) { int n = D * T; int p = XR.ncol(); if (D < 2) Rcpp::stop("D must be at least 2."); if (B < 1) Rcpp::stop("B must be positive."); // Map R input to Eigen Eigen::Map y(yR.begin(), n); Eigen::Map X(XR.begin(), n, p); Eigen::Map vardir(vardirR.begin(), n); // ---------- Fit the model to the original dataset ---------- VectorXd params(3); params[0] = 0.5 * vardir.mean(); // sig_u2 init params[1] = 0.5 * vardir.mean(); // sig_v2 init params[2] = 0.5; // rho init bool succ = fit_stfh_reml(params, y, X, vardir, D, T); if (!succ) Rcpp::stop("REML did not converge on the original data."); double sig_u2_hat = params[0]; double sig_v2_hat = params[1]; double rho_hat = params[2]; // ---------- Original EBLUP ---------- VectorXd theta_hat_orig(n); { MatrixXd J_T = MatrixXd::Constant(T, T, 1.0); MatrixXd C = ar1_corr(rho_hat, T); double inv_den = 1.0 / (1.0 - rho_hat * rho_hat); MatrixXd Omega = inv_den * C; MatrixXd InvSigmaSum = MatrixXd::Zero(p, p); VectorXd XtSigma_y = VectorXd::Zero(p); for (int d = 0; d < D; ++d) { VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2_hat * J_T + sig_v2_hat * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); MatrixXd XtVinv = Xd.transpose() * V_inv; InvSigmaSum += XtVinv * Xd; XtSigma_y += XtVinv * yd; } LLT lltK(InvSigmaSum); VectorXd beta_hat = lltK.solve(XtSigma_y); MatrixXd K_inv = lltK.solve(MatrixXd::Identity(p, p)); // Pre‑compute X * beta VectorXd Xbeta_full = X * beta_hat; for (int d = 0; d < D; ++d) { VectorXd Xbeta = Xbeta_full.segment(d * T, T); VectorXd yd = y.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2_hat * J_T + sig_v2_hat * Omega + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); VectorXd resid = yd - Xd * beta_hat; VectorXd Vinv_resid = V_inv * resid; // EBLUP components double u_hat = sig_u2_hat * VectorXd::Ones(T).dot(Vinv_resid); VectorXd v_hat = sig_v2_hat * Omega * Vinv_resid; theta_hat_orig.segment(d * T, T) = Xbeta + u_hat * VectorXd::Ones(T) + v_hat; } } // ---------- Generate all bootstrap random draws using R's RNG ---------- // This ensures full reproducibility Rcpp::NumericMatrix u_boot(D, B); Rcpp::NumericMatrix v_init_boot(D, B); Rcpp::NumericMatrix eta_boot(D * (T - 1), B); Rcpp::NumericMatrix e_boot(n, B); double sd_u = std::sqrt(sig_u2_hat); double sd_v1 = std::sqrt(sig_v2_hat / (1.0 - rho_hat * rho_hat)); double sd_v = std::sqrt(sig_v2_hat); for (int b = 0; b < B; ++b) { for (int d = 0; d < D; ++d) { u_boot(d, b) = R::rnorm(0.0, sd_u); v_init_boot(d, b) = R::rnorm(0.0, sd_v1); } for (int i = 0; i < D * (T - 1); ++i) { eta_boot(i, b) = R::rnorm(0.0, sd_v); } for (int i = 0; i < n; ++i) { e_boot(i, b) = R::rnorm(0.0, std::sqrt(vardirR[i])); } } // Pre‑compute X * beta_hat VectorXd Xbeta_full = X * lltK.solve(..)? No, we have beta_hat from original. // Actually compute Xbeta_full once MatrixXd InvSigmaSum_orig = /* recompute or store from above */; LLT lltK_orig(InvSigmaSum_orig); VectorXd beta_hat = lltK_orig.solve(XtSigma_y_orig); VectorXd Xbeta_full = X * beta_hat; // ---------- Bootstrap ---------- VectorXd mse_sum = VectorXd::Zero(n); VectorXd reblup_sum = VectorXd::Zero(n); VectorXd bias_sum = VectorXd::Zero(n); int n_failed = 0; #pragma omp parallel reduction(+ : n_failed) { VectorXd mse_local = VectorXd::Zero(n); VectorXd reblup_local = VectorXd::Zero(n); VectorXd bias_local = VectorXd::Zero(n); #pragma omp for for (int b = 0; b < B; ++b) { VectorXd y_b(n); VectorXd theta_b(n); // Construct bootstrap population and sample for (int d = 0; d < D; ++d) { double u_b = u_boot(d, b); VectorXd v_b(T); v_b[0] = v_init_boot(d, b); for (int t = 1; t < T; ++t) { int idx = d * (T - 1) + (t - 1); v_b[t] = rho_hat * v_b[t - 1] + eta_boot(idx, b); } VectorXd theta_linear = Xbeta_full.segment(d * T, T); VectorXd theta_d = theta_linear + u_b + v_b; theta_b.segment(d * T, T) = theta_d; for (int t = 0; t < T; ++t) { y_b[d * T + t] = theta_d[t] + e_boot(d * T + t, b); } } // Fit REML to bootstrap data VectorXd params_b(3); params_b[0] = sig_u2_hat; params_b[1] = sig_v2_hat; params_b[2] = rho_hat; bool success = fit_stfh_reml(params_b, y_b, X, vardir, D, T, 50); if (!success) { ++n_failed; continue; } // Compute bootstrap EBLUP (same structure as original) double sig_u2_b = params_b[0]; double sig_v2_b = params_b[1]; double rho_b = params_b[2]; MatrixXd C_b = ar1_corr(rho_b, T); double inv_den_b = 1.0 / (1.0 - rho_b * rho_b); MatrixXd Omega_b = inv_den_b * C_b; MatrixXd J_T = MatrixXd::Constant(T, T, 1.0); MatrixXd InvSigmaSum_b = MatrixXd::Zero(p, p); VectorXd XtSigma_y_b = VectorXd::Zero(p); for (int d = 0; d < D; ++d) { VectorXd yd_b = y_b.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); MatrixXd V = sig_u2_b * J_T + sig_v2_b * Omega_b + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); MatrixXd XtVinv = Xd.transpose() * V_inv; InvSigmaSum_b += XtVinv * Xd; XtSigma_y_b += XtVinv * yd_b; } LLT lltK_b(InvSigmaSum_b); VectorXd beta_b = lltK_b.solve(XtSigma_y_b); VectorXd Xbeta_b = X * beta_b; VectorXd theta_b_hat(n); for (int d = 0; d < D; ++d) { VectorXd yd_b = y_b.segment(d * T, T); MatrixXd Xd = X.middleRows(d * T, T); VectorXd psi = vardir.segment(d * T, T); VectorXd Xbeta_b_d = Xbeta_b.segment(d * T, T); MatrixXd V = sig_u2_b * J_T + sig_v2_b * Omega_b + psi.asDiagonal().toDenseMatrix(); LLT lltV(V); MatrixXd V_inv = lltV.solve(MatrixXd::Identity(T, T)); VectorXd resid = yd_b - Xd * beta_b; VectorXd Vinv_resid = V_inv * resid; double u_hat_b = sig_u2_b * VectorXd::Ones(T).dot(Vinv_resid); VectorXd v_hat_b = sig_v2_b * Omega_b * Vinv_resid; theta_b_hat.segment(d * T, T) = Xbeta_b_d + u_hat_b * VectorXd::Ones(T) + v_hat_b; } VectorXd diff = theta_b_hat - theta_b; for (int i = 0; i < n; ++i) { mse_local[i] += diff[i] * diff[i]; bias_local[i] += diff[i]; } reblup_local += theta_b_hat; } // end omp for #pragma omp critical { mse_sum += mse_local; reblup_sum += reblup_local; bias_sum += bias_local; } } // end omp parallel int B_eff = B - n_failed; VectorXd mse = mse_sum / B_eff; VectorXd reblup = reblup_sum / B_eff; VectorXd bias = bias_sum / B_eff; // ---------- Prepare output ---------- return Rcpp::List::create( Rcpp::Named("mse") = Rcpp::wrap(mse), Rcpp::Named("est") = Rcpp::List::create( Rcpp::Named("est") = Rcpp::wrap(theta_hat_orig), Rcpp::Named("reblup") = Rcpp::wrap(reblup), Rcpp::Named("bias") = Rcpp::wrap(bias) ) ); } ```

Rcpp code test