## 2 Materials and Methods
### Generalized Additive Model (GAM) Implementation
High \(\lambda\) (High Smoothness): When \(\lambda\) is large, the penalty for roughness dominates, forcing g to become smoother. At the extreme, \(s(. )\) reduces to a simple linear function (or a straight line). This results in low variance but high bias [40]. In this study, the optimal values of \(\lambda\) for s(x, y) and s(t) were chosen through leave-one-out cross-validation (LOOCV). This approach evaluates the performance of the model by minimizing the cross-validated residual sum of squares (RSS):  
\[\text{RSS}_{\text{cv}}(\lambda)=\sum\limits_{i=1}^{n}\left(\frac{y_{i}-g_{ \lambda}(x_{i})}{1-S_{\lambda,i}}\right)^{2} \tag{3}\]  
where \(S_{\lambda}\) is the smoothing matrix associated with the fit and \(S_{\lambda,i}\) represents its diagonal elements. LOOCV is computationally efficient for smoothing splines, as it avoids the need to refit the model repeatedly by leveraging the properties of \(S_{\lambda}\). A crucial byproduct of \(\lambda\) is its control over the effective degrees of freedom (EDF), which measure the complexity of the fitted function. For smoothing splines, the EDF is computed as:  
\[\text{EDF}=\sum\limits_{i=1}^{n}S_{\lambda,i} \tag{4}\]  
where a higher EDF indicates greater flexibility in the model. As \(\lambda\) increases, the EDF decreases, reducing the complexity of the fit. This study explicitly reports the EDF for each smoothing function to ensure interpretability and transparency in the modeling process [40]. For the spatial smoothing function \(s(x,\ y)\), the penalty term \(\lambda_{LatLon}\) controls smoothness along latitude and longitude values to ensure that the model captures key spatial patterns, such as high-demand regions for shelters, while avoiding overfitting to noise. A higher \(\lambda_{LatLon}\) results in smoother spatial surfaces, reducing the risk of overfitting to local noise, while a lower value allows the model to capture finer spatial details, such as localized demand hotspots. Similarly, for the temporal smoothing function \(s(t)\), the penalty term \(\lambda_{Day}\) controls the model's sensitivity to daily trends in vehicle counts. This parameter ensures a balance between identifying meaningful temporal patterns (e.g., evacuation peaks) and avoiding overfitting to random daily fluctuations. The smoothing parameters (\(\lambda\)) were optimized using Leave-One-Out Cross-Validation (LOOCV) approximated by the efficient Generalized Cross-Validation (GCV) criterion, as implemented in the mgcv package (version 1.9-1) in R [41]. The optimization focuses on minimizing the Penalized Residual Sum of Squares (PRSS) to balance model fit and smoothness. To ensure computational efficiency, the model uses Penalized Iteratively Reweighted Least Squares (P-IRLS) for convergence, which allows efficient optimization even with multiple smooth terms. Smoothing parameters for spatial (\(\lambda_{latlon}\)) and temporal (\(\lambda_{day}\)) effects were tuned to prevent overfitting while capturing key spatial-temporal patterns in evacuation behavior. _Prediction of demand in unobserved facilities:_ Once Equation (1) is estimated from data using statistical methods, the estimated coefficients, \(\beta_{0}\), \(\beta_{1}\),..., \(\beta_{7}\), \(\gamma\), penalty terms \(\lambda_{LatLon}\), \(\lambda_{Day}\) and smooth functions \(s(x,\ y)\) and \(s(t)\) will be available. Using these estimated quantities, the vehicle counts for a given spatial coordinate \((x,\ y)\) and day \((t)\), and the basis functions will constitute our predictions. _Model Evaluation:_ The performance of the model was evaluated using summary statistics, namely the Pearson's Correlation Coefficient (5), Root Mean Squared Error (6),Mean Absolute Error (7), and Mean Absolute Percentage Error (8), and visualizations, including actual vs. predicted train and test datasets plots.