--- title: "Demonstrating functional mediation" author: "John J. Dziak" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Demonstrating functional mediation} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` # Introduction The ``funmediation`` package fits a functional mediation model to a dataset consisting of intensive longitudinal data for a sample of individuals. For each individual $i$, the model assumes a time-invariant (i.e., non-time-varying) randomized treatment or exposure $X_i$, a distal outcome $Y_i$ observed at one time point (e.g., end-of-study), and a time-varying mediator $M_i(t)$ observed repeatedly during the interval after $X_i$ is observed and before $Y_i$ is observed. $X_i$ is treated as a number; that means it can either be dichotomous (and coded as 0 for no and 1 for yes) or continuous (and entered as the number). If $X_i$ is categorical with more than two categories, it needs to be entered as multiple binary dummy codes. The research question is whether the effect of $X_i$ on $Y_i$ is mediated by the process $M_i(t)$. A similar model was first proposed by Lindquist (2012). $M_i(t)$ can be continuous or binary, and $Y_i$ can also be either continuous or binary. The treatment variable $X_i$ is assumed to be binary. # Notation The mediation model for ``funmediation`` is fit in stages. First, the effect of $X_i$ on $M_i$ is fit behind the scenes as a time-varying effects (longitudinal varying coefficients) model (TVEM; see Hastie and Tibshirani, 1993; Tan et al., 2012). The assumed mean model is $E(M_i(t)) = \alpha_0(t) + \alpha_X(t) X_i$ for continuous $M_i(t)$, or $\mathrm{logit}^{-1}(E(M_i(t))) = \alpha_0(t) + \alpha_X(t) X_i$ for binary $M_i(t)$. This fits the marginal (population-averaged) mean: that is, without random effects but with a sandwich covariance estimate to handle within-subject correlation, as in working-independence generalized estimating equations. This part of the model is fit using the ``tvem`` function in the ``tvem`` R package. Next, the effect of $X_i$ and $M_i(t)$ on $Y_i$ is modeled as a scalar-on-function functional regression (see Goldsmith et al., 2011). The assumed mean model is either $E(Y_i) = \beta_0 + \beta_{X}X_i+ \int\beta_M(t)M_i(t)dt$ for continuous $Y_i$, or $\mathrm{logit}^{-1}(E(Y_i)) = \beta_0 + \beta_{X}X_i+ \int\beta_M(t)M_i(t)dt$ for binary $Y_i$. This part of the model is fit using the ``pfr`` function in the ``refund`` R package (see Goldsmith et al., 2011). Both the ``tvem`` and ``refund`` packages use the ``mgcv`` package (see Wood, 2017) for back-end calculations. In this mediation model, $\alpha_X(t)$ and $\beta_M(t)$ both must be nonzero in order for an indirect effect (i.e., mediation) to exist. These two functions arise conceptually from two different kinds of functional regression; one represents a model with a concurrent function as a response (of which TVEM is a special case), while the other represents a model with a distal scalar response (see Ramsay and Silverman, 2005). However, intuitively the size of the indirect effect has to do with both of them combined. Here we operationally define the indirect effect as the integral $\int\alpha_X(t)\beta_M(t)dt$. Because these functions are actually only estimated on a grid of points, the integral is approximated as a weighted average of the cross-products of the estimates. We obtain a bootstrap confidence interval for this quantity using the ``boot`` package. # Including Covariates Covariates can be included in predicting $M_i(t)$ and $Y_i$. For example, suppose there is a covariate $Z_i$ which has a time-varying relationship to $M_i(t)$. The TVEM model for the mediator can be expanded to $E(M_i(t)) = \alpha_0(t) + \alpha_X(t) X_i + \alpha_Z(t) Z_i$ for continuous $M_i(t)$, or $\mathrm{logit}^{-1}(E(M_i(t))) = \alpha_0(t) + \alpha_X(t) X_i+ \alpha_Z(t) Z_i$ for binary $M_i(t)$. It is possible for $Z_i$ to have a time-varying effect even if the values of $Z_i$ do not vary over time. That would mean that the correlation between the observation-level $M_i(t)$ and the subject-level $Z_i$ is different for different values of $t$. It is currently not *possible* to use a time-varying covariate in this package. But it is possible to assume that the relationship between $Z_i$ and $M_i(t)$ does not depend on $t$, whether or not $Z_i$ depends on $t$; in this case $\alpha_Z(t)$ can simply be written as $\alpha_Z$. The package allows both time-varying-effects and time-invariant-effects of covariates to be specified in predicting the mediator, using the ``tve_covariates_on_mediator`` and ``tie_covariates_on_mediator`` arguments, respectively. Alternatively, suppose that there is a subject-level covariate, $S_i$, which predicts $Y_i$. This can be added to the functional regression model by specifying $E(Y_i) = \beta_0 + \beta_{X}X_i+ \beta_S S_i + \int\beta_M(t)M_i(t)dt$ for continuous $Y_i$, or $\mathrm{logit}^{-1}(E(Y_i)) = \beta_0 + \beta_{X}X_i+ \beta_S S_i + \int\beta_M(t)M_i(t)dt$ for binary $Y_i$. Such a covariate can be included using the ``covariates_on_outcome``. Currently, the package assumes a subject-level $Y_i$ and does not support multiple functional coefficients in predicting the outcome; that is, the covariates in ``covariates_on_outcome`` cannot have time-varying values or effects. # Example with Continuous Outcomes The following example shows how to simulate example data and then analyze it using the ``funmediation`` package. ## Getting ready to run the example Before running the examples, first install and load the `funmediation` package. A `.zip` or `.tar.gz` file containing the package is available at [https://github.com/dziakj1/funmediation_development](https://github.com/dziakj1/funmediation_development), and it can then be used with the `install.packages()` function in R code, `Packages > Install Package(s)` from Local Files in the R graphical user interface, or `Tools > Install Packages` in the RStudio application, to install the package. We have also released `funmediation` on the CRAN archive, so that it can be installed using the command ```{r,eval=FALSE} install.packages("funmediation") ``` Another option is to install the package by code from the GitHub repository as follows: ```{r,eval=FALSE} install.packages("devtools") library(devtools) install_github("dziakj1/funmediation") ``` Of course, if you are viewing this guide from within R using the `vignette()` function, then the package is already installed. The next step is to load the required packages. ```{r} library(tvem) library(refund) library(boot) library(funmediation) ``` We then set a seed for replicability and simulate data. ```{r} set.seed(123) simulation1 <- simulate_funmediation_example(nsub=500) ``` The ``simulation1`` object will contain not only the simulated dataset itself, but the true values of the simulated parameters, including the indirect effect. ```{r} str(simulation1) ``` We need the simulated dataset as a data.frame object. ```{r} the_data <- simulation1$dataset ``` We can use the ``head`` and ``summary`` function in R, in order to look at the first few lines of the data and univariate descriptive statistics. ```{r} print(head(the_data)) summary(the_data) ``` It would be reasonable to do other descriptive analyses but in order to demonstrate the function we proceed immediately to the main analysis. # Running a Functional Mediation Model Now we call the functional mediation function. Only 10 bootstrap samples are used below for this quick illustration. At least a few hundred bootstraps are recommended in practice in order to increase precision and power. ```{r} model1 <- funmediation(data=the_data, treatment=X, mediator=M, outcome=Y, id=subject_id, time=t, nboot=10) ``` The warning message, ``extreme order statistics used as endpoints,'' occurs because too few bootstrap samples are used for the bootstrap confidence intervals to be valid; however, we ignore this here in order for the vignette example to run quickly. The `print` function gives an initial overview of the results. ```{r} print(model1) ``` The first part of the output above summarizes the estimated indirect effect and its bootstrap confidence interval. The estimate for the indirect effect $\int \alpha_X(t)\beta_M(t)$ is given as about -0248. Choosing the wider confidence interval in order to be conservative, a 95% confidence interval extends from about -0.35 to -0.14. (This interval would not actually be reliable in practice because of the low number of bootstrap samples, which could be remedied by choosing a higher ``nboot``). Thus, there appears to be a significant negative indirect effect of $X$ on $Y$ mediated through $M(t)$. The second part of the printed output summarizes the TVEM used to predict $M(t)$ from $X$. It only provides summary information on the model that was fit. The time-varying coefficients are not actually printed, because they are functions rather than single numbers. They can be plotted using the ``plot`` function, as described later. They are also stored in the ``model1`` output object, as described later. The third part summarizes the scalar-on-function functional regression model used to predict $Y$ from $X$ and $M(t)$. This model involves scalar (non-time-varying) coefficients for the intercept and the direct effect of $X$, and a functional (time-varying) coefficient for the effect of $M(t)$. The scalar coefficients are printed, but as before, the functional coefficient needs to be plotted using the ``plot`` function. The fourth section of the printed output simply presents an estimate of the total effect of $X$ on $Y$ ignoring $M(t)$, obtained using the ``glm`` function. $X$ is shown to have a total negative effect on $Y$, with an estimated coefficient of ``-0.30617`` and $p$-value of ``0.00175``. This suggests that the $X=1$ group has statistically significantly lower average $Y$ than the $X=0$ group, irrespective of $M$. The `plot` function plots the results. Three kinds of plots are available, represented by the options `tvem`, `pfr`, and `coef`. The `tvem` plot shows the estimates for $\alpha_0(t)$ and $\alpha_1(t)$. ```{r} plot(model1, what_plot="tvem") ``` The intercept plot appears to be significantly above zero and significantly increasing. It can be interpreted as the mean $M(t)$ for the control ($X=0$) group because $E(M(t)|X=0)=\alpha_0(t) + \alpha_X(t) \times 0=\alpha_0(t)$. The treatment effect plot is generally nonzero and decreasing. It can be interpreted as $E(M(t)|X=1)-E(M(t)|X=0)$, a time-specific treatment effect on the mediator, and it suggests that as time goes on the $X=1$ group tends to become lower in average $M(t)$ than the $X=0$ group. The `pfr` plot shows the estimate for $\beta_M(t)$. ```{r} plot(model1, what_plot="pfr") ``` The function tends to be nonzero at least for higher values of time, and generally seems to be steadily increasing. One possible interpretation is that the mediator becomes more important at later times. Alternatively, perhaps individuals who increase more (or decrease less) on average over time on the mediator tend to have a higher value on the outcome (see Dziak et al., 2019). The `pfrgam` plot does the same thing as the `pfr` plot, but with a slightly different implementation based on the `plot.gam` method from the `mgcv` library. ```{r} plot(model1, what_plot="pfrgam") ``` Last, the `coef` plot summarizes all of the most important coefficients estimated in the model. The upper left and upper right panes show the estimates of $\alpha_0(t)$ and $\alpha_X(t)$ from the model of the effect of $X$ on $M(t)$. The lower left shows the estimate of $\beta_M(t)$ from the model of the effect of $M(t)$ on $Y$ given $X$. In the lower right panel, the estimate of the indirect effect is printed (not plotted, because it is a single number). ```{r} plot(model1, what_plot="coef") ``` Because we have already looked at the relevant plots before using the ``coef`` option, there is nothing new to interpret in the ``coef`` plot; it is just intended as a convenient summary. ## Exploring Information in the Output Objects The output object created, called ``model1``, contains useful information both on the model fit to the original data and on the bootstrapping results. * ``model1$original_results$time_grid`` is the grid of time points on which the coefficient functions were estimated. * ``model1$original_results$alpha_int_estimate``, ``model1$original_results$alpha_int_se``, ``model1$original_results$alpha_X_estimate``, and ``model1$original_results$alpha_X_se`` are the functional coefficients and pointwise standard errors for the intercept and treatment in the TVEM, at each point of the time grid. * ``model1$original_results$beta_int_estimate``, ``model1$original_results$beta_int_se``, ``model1$original_results$beta_X_estimate`` and ``model1$original_results$beta_X_se`` are the estimates and standard errors for the scalar coefficients representing the intercept and the direct treatment effect in the scalar-outcome functional linear model. * ``model1$original_results$beta_M_estimate`` and ``model1$original_results$beta_M_se`` are the estimates and pointwise standard errors representing the effect of the time-varying mediator at each point of the time grid. * ``beta_M_pvalue`` is the $p$-value for an overall test of whether the mediator is significantly associated with the outcome, i.e., whether $\beta_M(t)$ can be shown not to be all zeroes; see the documentation for the ``pfr`` function in the ``refund`` package for more information. * ``model1$original_results$tau_int_estimate``, together with ``$tau_int_se``, ``$tau_X_estimate``, and ``$tau_X_se``, represent the estimates of the intercept and coefficient in the model for the total effect of $X$ on $Y$. Thus, ``model1$original_results$tau_X_estimate`` is the estimated total effect and ``model1$original_results$tau_X_pvalue`` is the $p$-value for testing whether this effect is zero. * ``model1$original_results$indirect_effect_estimate`` gives the estimated value of $\int \alpha_X(t)\beta_M(t)$ obtained without bootstrapping. * ``model1$original_results$tvem_XM_details``, ``$funreg_MY_details``, and ``$total_effect_details`` give the fitted model objects from ``tvem``, ``pfr`` and ``glm``. They might be useful for plots or diagnostics but might contain more technical detail than some users are interested in. ## Including a Covariate It is possible to run the model with one or more additional covariates. To demonstrate this, the ``simulate_mediation_example`` has an option ``make_covariate_S`` to create an extra subject-level covariate. The data is generated in such a way that $S$ is not actually related to $X$, $M$ or $Y$, but the analyst is assumed not to know this -- that is, the simulated true value of its coefficient is zero. ```{r} simulation_with_covariate <- simulate_funmediation_example(nsub=500, make_covariate_S=TRUE); data_with_covariate <- simulation_with_covariate$dataset; ``` To include the covariate $S_i$ with an assumption of a possibly time-varying effect $\alpha_S(t)$ on $M(t)$ (recall that the covariate can have a time-varying effect even though it may not have time-varying values), use the `covariates_on_outcome` and `tve_covariates_on_mediator` arguments: ```{r} model_with_tve_covariate <- funmediation(data=data_with_covariate, treatment=X, mediator=M, outcome=Y, tve_covariates_on_mediator = ~S, covariates_on_outcome = ~S, id=subject_id, time=t, nboot=10) ``` Notice that the tilde (~) character needs to be included for the covariates. This is done so that R will allow multiple covariates to be listed in the style of the right side of a formula, e.g., ``~S1 + S2``. The covariates in ``tve_covariates_on_mediator`` and ``covariates_on_outcome`` need to be subject-level and should not have time-varying values; the function does not currently support models with time-varying values other than the mediator. However, the relationship of the covariates to the mediator can be specified as either time-varying or time-invariant. The relationship of the covariate to the outcome, in contrast, can only be specified as time- invariant in the current version of the package. That is, $M(t)$ is the only time-varying predictor in the functional regression predicting the outcome. The results of the model fit can be viewed using the `print` and `plot` functions: ```{r} print(model_with_tve_covariate) plot(model_with_tve_covariate, what_plot="tvem") ``` Notice that we now get a time-varying effects plot which includes a coefficient for the effect of $S$, in addition to the coefficients for the effect of $X$. In this example $\alpha_S(t)$ is not significantly nonzero at any time and does not change over time, so there is no evidence that covariate $S$ is important in predicting the outcome $Y$. As before, the $X=1$ group has lower values on the mediator than the $X=0$ group (i.e., their contrast is negative) at least at later time points. Alternatively, by specifying ``tie_covariates_on_mediator`` instead of ``tve_covariates_on_mediator,`` the covariate can be assumed to have a time-invariant effect rather than a time-varying effect on $M$. ```{r} model_with_tie_covariate <- funmediation(data=data_with_covariate, treatment=X, mediator=M, outcome=Y, tie_covariates_on_mediator = ~S, covariates_on_outcome = ~S, id=subject_id, time=t, nboot=10) print(model_with_tie_covariate) plot(model_with_tie_covariate, what_plot="pfr") ``` The functional effect plot here is very similar to what it was before. Now there are scalar estimates for $\alpha_S=-0.03008$ and $\beta_S=-0.00164$ but no plot for $\alpha_S(t)$ because $\alpha_S$ is now assumed not to vary over time. Again, $S$ is not statistically significant either as a predictor of the mediator or of the outcome, because the coefficients are not large relative to their standard errors. ## Using a Binary Mediator The mediator may be binary rather than continuous, in which case the the TVEM uses the logistic link function. For purposes of this vignette demonstration, we dichotomize $M$ in the simulated dataset and use that as our binary mediator (this usually would not be advisable for real data because it involves loss of information). ```{r} data_with_covariate$binary_M <- 1*(data_with_covariate$M > mean(data_with_covariate$M)) head(data_with_covariate) ``` Now we call the `funmediation` function and include the `binary_mediator=TRUE` argument. ```{r} model_with_binary_M <- funmediation(data=data_with_covariate, treatment=X, mediator=binary_M, outcome=Y, id=subject_id, time=t, binary_mediator=TRUE, nboot=10) ``` ```{r} print(model_with_binary_M) plot(model_with_binary_M, what_plot="coef") ``` The plots are similar to those shown in the previous example with numerical mediator, except that the TVEM coefficient functions for the intercept and treatment effects on the outcome are interpreted on the logit scale as in logistic regression. Note that if the ``binary_mediator=TRUE`` argument is not specified, the TVEM will be fit using an identity link (hence a linear model at any given time $t$) as for a normally distributed outcome. The function would not detect that ``binary_M`` consists of zeroes and ones and no other numbers. The identity link is probably not the best approach but might be interesting as a comparison. If the mediator is binary, it is probably better to use ``binary_mediator=TRUE`` in order to specify that the $\alpha$ coefficients are to be interpreted on the logit (logistic regression) scale. A future version of this package might allow log link functions as an alternative to logit link functions, but this is currently not available. # Example with a Binary Outcome In addition to, or instead of, a binary mediator, the outcome may also binary. In this case, the scalar-on-function part of the model (the $\beta_0$, $\beta_X$ and $\beta_M(t)$ functions) are fit using a logistic link function. To demonstrate this, we first simulate a dataset with a binary outcome, using the ``simulate_binary_Y`` argument in ``simulate_funmediation_example`` function. We set this argument to TRUE, as it is FALSE by default. (Note that we do not currently have a corresponding ``simulate_binary_M`` argument in this function, as simulating realistic binary data trajectories is somewhat more difficult.) ```{r} simulation_with_binary_Y <- simulate_funmediation_example(simulate_binary_Y=TRUE, nsub=500) data_with_binary_Y <- simulation_with_binary_Y$dataset head(data_with_binary_Y) ``` Now we fit the model, using the new dataset with binary $Y$, and using the ``binary_outcome=TRUE`` argument. By default, ``binary_outcome`` is FALSE. The ``binary_outcome=TRUE`` argument specifies that the $\beta$ coefficients are to be interpreted as logistic regression coefficients, similar to the ``binary_mediator=TRUE`` argument which does the same for the $\alpha$ coefficients. ```{r} model_with_binary_Y <- funmediation(data=data_with_binary_Y, treatment=X, mediator=M, outcome=Y, id=subject_id, time=t, binary_outcome=TRUE, nboot=10) ``` ```{r} print(model_with_binary_Y) plot(model_with_binary_Y, what_plot="coef") ``` Note that the relationship of the mediator function to the outcome is now interpreted as a functional logistic regression coefficient rather than a functional linear regression coefficient (see Dziak et al., 2019; Goldsmith, Crainiceanu, Caffo & Rice, 2011; Mousavi and S\o rensen, 2018 for descriptions of functional logistic regression). # More than Two Treatment Groups In the preceding examples, the treatment was assumed to be dichotomous and dummy-coded (0 for a control group and 1 for a treated group). Thus, the coefficients were interpreted as effects of increasing the binary variable $X_i$ from 0 to 1. However, sometimes in practice there might be more than two groups in the experiment; for example, there might be two different interventions each being compared to the same control group. In this case, $X_i$ should not be coded as, say, 1, 2, and 3, for the different groups, because in that case the `funmediation` function would misinterpret the levels as if they were actually numbers. To handle this, dummy-code multiple predictor variables so that, e.g., `X1` is 1 for treatment group 1 and 0 otherwise, and `X2` is 1 for treatment group 2 and 0 otherwise. In this case, the effect of `X1` represents the contrast between group 1 and the remaining (reference) group 3, and the effect of `X2` represents the contrast between group 2 and the reference group. You could optionally choose a different reference group (e.g., have a dummy code for group 2 and a dummy code for group 3, leaving 1 as the reference level, instead of a dummy code for groups 1 and 2). Except when using data simulated by the package, you should code the predictor variables yourself, keeping in mind which contrasts are most interesting to you. Regardless, one of the treatment levels, perhaps representing the control or reference group, does not get its own predictor variable but is represented indirectly by having zeroes on all the others. The following code will simulate data from a study with three treatment groups, hence two binary treatment variables. ```{r} set.seed(12345) nboot <- 99 answers <- NULL f1 <- function(t) {return(-(t/2)^.5)} f2 <- function(t) {return(sqrt(t))} f3 <- function(t) {return(sin(2*3.141593*t))} the_simulation <- simulate_funmediation_example(nlevels=3, alpha_X = list(f1,f2), beta_M = f3, beta_X = c(.2,.3)) ``` In the code above, notice that `alpha_X`, the time-varying effect of treatment on the mediator, has to be specified as a list of two functions, corresponding to the two dummy-coded dimensions of treatment. Similarly, `beta_X`, the direct effect of treatment on outcome, is a vector of two numbers. The first few lines of the resulting dataset look like this: ```{r} head(the_simulation$dataset) ``` The following code will analyze this data: ```{r} the_data <- the_simulation$dataset ans_funmed_3groups <- funmediation(data=the_data, treatment=~X1+X2, mediator=M, outcome=Y, tvem_num_knots=5, id=subject_id, time=t, nboot=10) ``` The `print` function works as before to print the model details. ```{r} print(ans_funmed_3groups) ``` The `pfr` plot is essentially the same as before, since there is only one mediator. ```{r} plot(ans_funmed_3groups, use_panes=TRUE, what_plot = "pfr") ``` In this example, the mediator seems to be positively related to the outcome at earlier times, but weakly negatively related at later times, although substantive interpretation of such functions requires care. The `tvem` plot now shows the time-varying effects for both predictors. ```{r} plot(ans_funmed_3groups, use_panes=TRUE, what_plot = "tvem") ``` The coefficient for the intercept shows that at least for the reference group (group 3 in this example) the expected value of the response tends to increase over time. The coefficient for the first treatment variable (contrasting group 1 vs. group 3) shows a lower value on the mediator for group 1 relative to group 3 (a negative contrast) at least at the later time points. The coefficient for the second treatment variable suggests that group 2 has a higher mean than group 3, especially at later time points. The `coefs` plot also includes the extra time-varying effect. The indirect effects are not printed in the lower right pane anymore, because that pane is needed for the new time-varying effect. However, it can be viewed in the output from the `print` statement as usual. ```{r} plot(ans_funmed_3groups, use_panes=TRUE, what_plot = "coefs") ``` Currently the package can still only handle a single mediator variable. It doesn't currently allow, for example, different mediator variables for different treatment variables. # References * Dziak, J. J., Coffman, D. L., Reimherr, M., Petrovich, J., Li, R., Shiffman, S., Shiyko, M. P. (2019). Scalar-on-function regression for predicting distal outcomes from intensively gathered longitudinal data: Interpretability for applied scientists. Statistical Surveys, 13: 150-180. * Goldsmith, J., Bobb, J., Crainiceanu, C., Caffo, B., and Reich, D. (2011). Penalized functional regression. Journal of Computational and Graphical Statistics, 20(4), 830-851. * Goldsmith, J., Crainiceanu, C. M., Caffo, B. S., & Reich, D. S. (2011). Penalized functional regression analysis of white-matter tract profiles in multiple sclerosis. Neuroimage, 57(2): 431–439. * Hastie, T., & Tibshirani, R. (1993). Varying-coefficient models. Journal of the Royal Statistical Socety, B, 55:757-796. * Lindquist, M. A. (2012). Functional Causal Mediation Analysis With an Application to Brain Connectivity. Journal of the American Statistical Association, 107: 1297-1309. * Ramsay, J. O., & Silverman, B. W. (2005). Functional data analysis (2nd ed.). New York: Springer. * Mousavi, N., & S\o rensen, H. (2018). Functional logistic regression: a comparison of three methods. Journal of Statistical Computation and Simulation, 88:2, 250-268. * Tan, X., Shiyko, M. P., Li, R., Li, Y., & Deirker, L. (2012). A time-varying effect model for intensive longitudinal data. Psychological Methods, 17: 61-77. * Wood, S.N. (2017) Generalized Additive Models: an introduction with R (2nd edition), CRC. # Acknowledgements The ``pfr`` function was developed by Jonathan Gellar, Mathew W. McLean, Jeff Goldsmith, and Fabian Scheipl, the ``boot`` package was written by Angelo Canty and Brian Ripley, the ``mgcv`` package was developed by Simon N. Wood, and the ``refund`` package was developed and maintained by Julia Wrobel and others.