Abstract
We propose a novel modeling framework to study the effect of covariates of various types on the conditional distribution of the response. The methodology accommodates flexible model structure, allows for joint estimation of the quantiles at all levels, and provides a computationally efficient estimation algorithm. Extensive numerical investigation confirms good performance of the proposed method. The methodology is motivated by and applied to a lactating sow study, where the primary interest is to understand how the dynamic change of minute-by-minute temperature in the farrowing rooms within a day (functional covariate) is associated with low quantiles of feed intake of lactating sows, while accounting for other sow-specific information (vector covariate).
1. Introduction
Many modern applications routinely collect data on study participants comprising scalar responses and covariates of various types, vector, function, and image, and the main question of interest is to examine how the covariates affect the response. For example, in our motivating experimental study, the goal is to analyze how the minute-by-minute daily temperature and humidity of the farrowing rooms, where sows are placed after giving birth for nursing, affect their feed intake during a lactation period. The covariates consist of temperature profile, humidity, and sow age, where the response is a total amount of daily feed intake of sows. A popular approach in these cases is to use a nonparametric framework and assume that the mixed covariates solely affect the mean response; see Cardot et al. [1], James [2], Ramsay and Silverman [3, 4], Ferraty and Vieu [5, 6], Goldsmith et al. [7], McLean et al. [8], and others. However, for our application, while it is important to study the mean feed intake, animal scientists are often more concerned with the left tail of the feed intake distribution. This is because low feed intake of lactating sows could lead to many serious issues, including decrease in milk production and negative impact on the sows reproductive system; see, for reference, Quiniou and Noblet [9], Renaudeau and Noblet [10], St-Pierre et al. [11], among others. In this paper we focus on regression models that study the effects of covariates on the entire distribution of response. Our contribution is the development of a modeling framework that accommodates a comprehensive study of various types (vector and functional) of covariates on a scalar response.
Quantile regression models the effect of scalar/vector covariates beyond the mean response; it provides a more comprehensive study of the covariates on the response and has attracted great interest [12, 13]. For prespecified quantile levels, quantile regression models the conditional quantiles of the response as a function of the observed covariates; this approach has been extended more recently to ensure noncrossing of quantile functions [14]. Quantile regression has been also extended to handle functional covariates. Cardot et al. [15] discussed quantile regression models by employing a smoothing spline modeling based approach. Kato [16] considered the same problem and used a functional principal component (fPC) based approach. Both papers mainly discussed the case of having a single functional covariate and it is not clear how to extend them to the case where there are multiple functional covariates or mixed covariates (vector and functional).
More recently, Tang and Cheng [17], Lu et al. [18], and Yu et al. [19] studied quantile regression when the covariates are of mixed types and introduced the partial functional linear quantile regression modeling framework. The first two publications used fPC basis while the last one considered partial quantile regression (PQR) basis. These approaches all are suitable when the interest is studying the effect of covariate at a particular quantile level and do not handle the study of covariate effects at simultaneous quantile levels due to the well-known crossing-issue.
Ferraty et al. [20] and Chen and Müller [21] considered a different perspective and studied the effect of a functional predictor on the quantiles of the response by modeling the conditional distribution of the response directly. However their approach is limited to one functional predictor. In this paper we fill this gap and propose a unifying modeling framework and estimation technique that allows studying the effect of mixed type covariates (i.e., scalar, vector, and functional) on the conditional distribution of a scalar response in a computationally efficient manner.
Let denote the th conditional quantile of given a functional covariate , and let denote the conditional distribution of given . We model the conditional distribution using a generalized function-on-function regression framework, i.e., , where is an indicator function and is the logit link function, and we study the conditional quantiles by exploiting the relationship between and through for . The advantage and contribution of our proposed method mainly come from the following reasons: (1) our modeling approach is spline-based, and as a result it can easily accommodate smooth effects of scalar variables as well as of functional covariates and (2) our estimation approach is based on a single step function-on-function (or function-on-scalar) penalized regression, which enables efficient implementation by exploiting off-the-shelf software and leads to competitive computations.
The remainder of the paper is structured as follows. Section 2 discusses the details of the proposed method and Section 3 describes the estimation procedure and extensions. Section 4 performs a thorough simulation study evaluating the performance of the proposed method and its competitors. We apply the proposed method to analyze the sow data in Section 5. We conclude the paper with a discussion in Section 6.
2. Methodology
2.1. Statistical Framework
Let be index subjects, index repeated measurements, the number of subjects, and the number of observations for subject . Suppose we observe , where is the response, which is a scalar random variable, is a -dimensional vector of nuisance covariates, is a scalar covariate, and is a functional covariate, which is assumed to be square-integrable on a closed domain .
We propose the following model for the conditional distribution of given , , and :where denotes the conditional distribution function as before, is a known, monotone link function, namely, the logit link function defined as for arbitrary scalar , is an unknown and smooth functional intercept, is a -dimensional parameter capturing the linear additive effect of the covariate vector , is an unknown and smooth function, and is an unknown and smooth bivariate function. Here, the effect of the nuisance covariates is ; it is assumed to be constant over while the smooth intercept is -variant. The effect of is , which varies smoothly over ; quantifies -variant linear effects of the covariate . If the parameter function is zero then the covariate has no effect on the distribution of the response , which is equivalent to having no effect on any quantile level of . Similarly, it is easy to see that a null effect, say , is equivalent to the case that the functional covariate has no effect on any quantile level of the response. Chen and Müller [21] (CM, henceforth) considered a similar model; however, their approach is restrictive to a single functional covariate. We discuss the differences between their method and ours in Section 2.2
To explain our ideas, we consider the case that the functional covariates are observed without noise on a fine, regular, and common grid of sampling points, i.e., with for all . Bear in mind, this assumption is made for illustration only, and our framework can be extended to more general cases, including settings where the functional covariate is observed with noise and at irregular sampling points; see Section 3.3.
2.2. Modeling of the Covariate Effects
We model and by using prespecified, truncated univariate basis. Let and be two bases of dimensions and , respectively. and , where ’s and ’s are unknown basis coefficients. We represent using the tensor product of two univariate bases functions, and , where and are the bases dimensions; , where ’s are unknown basis coefficients. In practice, the integration term is approximated by Riemann integration but other numerical approximation scheme can be also used.
Define for . In practice for each in a fine grid, we view as a binary-valued random functional variable. It follows that model (1) can be written equivalently as a generalized function-on-function regression model through relating the “artificial” binary functional response and the mixed covariates . This model can be fitted by using, for example, the ideas of Scheipl et al. [22] which we briefly summarize next.
Model (1) can be represented as the following generalized additive model:For convenience, we use the notation , . We let , , , , , , and is a coefficient matrix. Then , where is the vectorization of . We let , , and . Now model (2) can be written asThe general idea is to set the bases dimensions , , , and to be sufficiently large to capture the complexity of the coefficient functions and control the smoothness of the estimator through some roughness penalties. This approach of using roughness penalties has been widely used; see, for example, Eilers and Marx [23]; Ruppert [24]; Wood [25, 26] among many others.
It is important to emphasize that, even in the case of a single functional covariate, our methodology differs from [21] in two directions: (1) our proposed method is based on modeling the unknown smooth coefficient functions using prespecified basis function expansion and using penalties to control their roughness. In contrast, CM uses data-driven basis and chooses the number of basis functions through the percentage of explained variance (PVE) of the functional predictors. This key difference allows our method to accommodate covariates of different types as well as nonlinear effects. (2) Our estimation approach is based on a single step penalized function-on-function regression while CM uses pointwise estimation based on functional principal component bases and thus requires fitting multiple generalized regressions. This nice feature leads to an computational advantage.
3. Estimation
3.1. Estimation via Penalized Log-Likelihood
Let be a set of equally spaced points in the range of the response variable, ’s. Conditioning on , we model as independently distributed Bernoulli variables with mean , where . The coefficients , , , and are estimated by minimizing the penalized log-likelihood criterion:where is the log-likelihood function of data , , , , and are smoothing parameters, which control the balance between the model fit and its complexity, and , , , and are all penalties.
There are several choices to define the penalty matrix in nonparametric regression; see Eilers and Marx [23] and Wood [26]. We use quadratic penalties which penalize the size of the curvature of the estimated coefficient functions. Let , where is of dimension with its element equal to . Similarly, , where is of dimension with its element equal to . As is a bivariate function, the choice of penalty implies penalizing the size of curvature in each direction, respectively: , where is of dimension with the element of equal to for some orthonormal spline bases. Similarly, , where is of dimension with the element of equal to .
Criterion (4) can be viewed as a penalized quasilikelihood (PQL) of the corresponding generalized linear mixed modelwhere is the generalized inverse matrix of ; and are defined similarly. Wood [26] discusses an alternative way to deal with the rank-deficient matrices in the context of restricted maximum likelihood (REML) estimation. Here we do not account for the dependence over ; see Scheipl et al. [22] for a general formulation. See also Ivanescu et al. [27] who use the mixed model representation of a similar regression model to (5), but with a Gaussian functional response. The smoothing parameters are estimated using REML.
3.2. Extension to Nonlinear Model
One advantage of the proposed framework is that it can be easily extended to allow for more flexible effects, i.e., extending the ideas to accommodate multiple covariates, scalar or functional, and varied types of effects. In particular, the smooth effect can be replaced by and by , where and are unknown bivariate and trivariate smooth functions, respectively; see Scheipl et al. [22] and Kim et al. [28]. These changes require little additional computational effort. The modeling and estimation follow roughly similar ideas as Scheipl et al. [22]. We consider the nonlinear model in the simulation study for the case of having a scalar covariate only, i.e., , and the corresponding results are presented in Section S1.1 of the Supplementary Materials. The results show excellent prediction performance compared to the competitive nonlinear quantile regression method, namely, constrained B-spline smoothing (COBS) [29].
3.3. Extension to Sparse and Noisy Functional Covariates
In practice the functional covariates are often observed at irregular times across the units and the measurements are possibly corrupted by noises. In such case, one needs to first smooth and denoise the trajectories before fitting. When the sampling design of the functional covariate is dense, the common approach is to smooth each trajectory using splines or local polynomial smoothing, as proposed in Ramsay and Silverman [3] and Zhang and Chen [30]. When the design is sparse, the smoothing can be done by pooling all the subjects and following the PACE method proposed in Yao et al. [31]. As recovering the trajectories has been extensively discussed in the literatures, we do not review the procedures here. Instead, we discuss some available computing resources that can be used to fit these methods. In our numerical study, we used fpca.sc function in the refund R package [32] for recovering the latent trajectories, irrespective of a sampling design (dense or sparse). Alternatively, one can use fpca.face [33] in refund for regular dense design and face.spare [34] in the R package face [35] for irregular sparse design. Once the latent trajectories are estimated, they can be used in the fitting criterion (4).
3.4. Estimation of Conditional Quantile
Let , , , and be parameter estimates in (5). It follows that the estimated distribution function can be obtained by plugging in the estimated coefficients. The th conditional quantile is estimated by inverting the estimated distribution, i.e., . The estimated distribution function is not a monotonic function yet. In practice we suggest to first apply a monotonization method as described in Section 3.5 and then estimate the conditional quantiles by inverting the resulting estimated distribution.
3.5. Monotonization and Implementation
While a conditional quantile function is nondecreasing, the resulting estimated quantiles may not be. Two approaches are widely used: one is to monotonize the estimated conditional distribution function and the other is to monotonize the estimated conditional quantile function. We choose the former as is readily available at dense grid points ’s. We use an isotonic regression model [36] for monotonization, which imposes an order restriction; this is done by using the R function isoreg. Other monotonization approaches include Chernozhukov et al. [37], which was employed in Kato [16].
Our approach is implemented by first creating an artificial binary response and then fitting a penalized function-on-function regression model and using the logit link function. Fitting models in (4) can be done by extending the ideas of Ivanescu et al. [27] for Gaussian functional response; the extension of the model to the non-Gaussian functional response has recently been studied and implemented by Scheipl et al. [22] as the pffr function in refund package [32].
4. Simulation Study
4.1. Simulation Setting
In this section we evaluate the empirical performance of the proposed method. We present results for the case when we have both functional and scalar covariates; additional results when there is only a single scalar or a single functional covariate are discussed in the Supplementary Materials, Section S1.
Suppose the observed data for the th subject are , , where , for . Let , where , for odd values of , for even values of , , , and . We assume three cases for generating response :(i)Gaussian: ; this corresponds to the quantile regression model , where is the distribution function of the standard normal;(ii)Mixture of Gaussians: + , where the true quantiles can be approximated numerically by using qnorMix function in the R package norMix;(iii)Gaussian with heterogeneous error: ; the true quantiles are given by = + .
Let , where for .
For each case, we use different combinations of signal to noise ratio (SNR), sample size, and sampling designs to generate simulated datasets. We define SNR as , and we consider five levels of noise: . Two levels of sample size are and . Two sampling designs are considered: (i) sparse design, where are randomly selected points from a set of equispaced grids in ; and (ii) dense design, where the sampling points are equispaced time points in .
The performance is evaluated on a test set of subjects, for which we have available, in terms of mean absolute error (MAE) for quantile levels , and :
4.2. Competing Methods
We denote the proposed method by Joint QR to emphasize the single step estimation approach. We compare our method with two alternative approaches: (1) a variant of our proposed approach using pointwise estimation, denoted by Pointwise QR. This approach consists of fitting multiple regression models with binomial link function as implemented by the penalized functional regression pfr, developed by Goldsmith et al. [7], of the refund package for generalized scalar responses; (2) a modified version of the CM method, denoted by Mod CM, that we developed to account for additional scalar covariates and which fits multiple generalized linear models with scalar covariates and fPC scores as predictors; (3) a linear quantile regression approach using the quantile loss function and the partial quantile regression bases for functional covariates, proposed by Yu et al. [19] and denoted by PQR. Notice that although the formulation of the first two methods implicitly accounts for a varying effect of the covariates on the response distribution, they do not ensure that this effect is smooth. The third approach can only estimate a specific quantile rather than the entire conditional distribution. Note that all the competing methods are monotonized for a fair comparison.
The R function pfr can incorporate both scalar/vector and functional predictors by adopting a mixed effects model framework. The functional covariates are presmoothed by fPC analysis [31]. Throughout the simulation study we fix PVE as for fPC analysis to determine the number of principal components and use REML to select the smoothing parameters for our proposed methods. Other basis settings are set to their default values. We use 100 equally distanced points between the minimum and maximum of the observed ’s to set the grid for the conditional distribution function.
4.3. Simulation Results
Tables 1 and 2 show the accuracy of the quantile prediction for the two cases (normal and mixture) when the functional covariate is observed sparsely and the sample size is (Table 1) and (Table 2). Table 3 presents the prediction accuracy for the case of heteroskedasticity with sparsely observed functional covariates. The results based on dense sampling design show similar patterns and thus are relegated to the supplement; see Section S1.3. The comparison of running times is presented in Table 4.
For the case when the response is Gaussian, Tables 1 and 2 suggest that the Joint QR typically outperforms its competitors especially for lower quantile levels ( and ). For very small noise level (SNR = 150), PQR performs the best, followed closely by the proposed Joint QR. The variant Pointwise QR, which has a poorer performance, is generally better than the modified CM approach. As expected, as the sample size increases (), all the accuracy results improve; the proposed Joint QR continues to yield most accurate quantiles for the low quantile levels. For mixture of Gaussians, the results are somewhat similar. The accuracy of the quantile estimators with the Pointwise QR improves greatly; in fact the Joint QR and Pointwise QR outperform the other approaches for quantile levels irrespective of the SNR. Finally, Table 3 shows that the results for Gaussian with heterogeneous error are close to those for the case of Gaussian. Again, the proposed method has competitive performance in terms of prediction accuracy.
Table 4 compares the three methods that involve estimating the conditional distribution in terms of the running time required for fitting. The times are reported based on a computer with a 2.3 GHz CPU and GB of RAM. Not surprisingly by fitting the model a single time, Joint QR is the fastest, in some cases being order of magnitude faster than the rest. Pointwise QR can be up to twice as fast as Mod CM.
For completeness, we also compare our proposed method to the appropriate competitive methods for the cases (1) when there is a single scalar covariate and (2) when there is a single functional covariate. In the Supplementary Materials, Section S1.1 discusses the former case and compares Joint QR and Pointwise QR with the linear quantile regression (LQR) [12], implemented by rq function in the R package quantreg, and the constrained B-splines nonparametric regression quantiles (COBS), implemented by the cobs function in the R package COBS [29], in an extensive simulation experiment that involves both linear quantile settings and nonlinear quantile settings. Overall the results show that the proposed methods have similar behavior as LQR; see Table S1. Furthermore we consider the proposed method and its variant with nonlinear modeling of the conditional distribution as discussed in Section 3.2, which we denote with Joint QR (NL) for joint fitting and Pointwise QR (NL) for pointwise fitting. Nonlinear versions of the proposed methods have an excellent MAE performance, which is comparable to or better than that of the COBS method.
Finally, Section S1.2 in the Supplementary Materials discusses the simulation study for the case of having a single functional covariate and compares the proposed methods with CM in terms of MAE as well as computation time; see results displayed in Tables S2 and S3. The results show that the proposed Joint QR is comparable to CM in terms of the prediction accuracy and has less computation time. In our simulation study we also consider the joint fitting of the model by treating the binary response as normal and use pffr [27] with Gaussian link, denoted by Joint QR (G).
5. Sow Data Application
Our motivating application is an experimental study carried out at a commercial farm in Oklahoma from July 21, 2013, to August 19, 2013 [38]. The study comprises 480 lactating sows of different parities (i.e., number of previous pregnancies, which serves as a surrogate for age and body weight) that were observed during their first 21 lactation days; their feed intake was recorded daily as the difference between the feed offer and the feed refusal. In addition the study contains information on the temperature and humidity of the farrowing rooms, each recorded at five minute intervals. The final dataset we used for the analysis consists of 475 sows after five sows with unreliable measurements were removed by the experimenters.
The experiment was conducted to gain better insights into the way that the ambient temperature and humidity of the farrowing room affect the feed intake of lactating sows. Previous studies seem to suggest a reduction in the sow’s feed intake due to heat stress: above 29°C sows decrease feed intake by 0.5 kg per additional degree in temperature [9]. Studying the effect of heat stress on lactating sows is a very important scientific question because of a couple of reasons. First, the reduction of feed intake of the lactating sows is associated with a decrease in both their bodyweight (BW) and milk production, as well as the weight gain of their litter [10, 39, 40]. Sows with poor feed intake during lactation continue the subsequent reproductive period with negative energy balance [41], which leads to preventing the onset of a new reproductive cycle. Second, heat stress reduces farrowing rate (number of sows that deliver a new litter) and number of piglets born [42]; the reduction in reproduction due to seasonality is estimated to cost 300 million dollars per year for the swine industry [11]. Economic losses are estimated to increase [43] because high temperatures are likely to occur more frequently due to global warming [44].
Our primary goal is to understand the thermal needs of the lactating sows for proper feeding behavior during the lactation time. We are interested in how the interplay between the temperature and humidity of the farrowing room affects the feed intake demeanor of lactating sows of different parities. We focus on three specific time points during the lactation period—beginning (lactation day 4), middle (day 11), and end (day 18)—and the analyses are done separately for each time point. We consider two types of responses that are meant to assess the feed intake behavior using the current and the previous lactation days. The first one quantifies the absolute change in the feed intake over two consecutive days and the second one quantifies the relative change and takes into account the usual sow’s feed intake. We define them formally after introducing some notation.
Let be the th measurement of the feed intake observed for the th sow and denote by the lactation day when is measured; here . Most sows are observed every day within the first 21 lactation days and thus have . First define the absolute change in the feed intake between two consecutive days as for that satisfies . For instance, means there was no change in feed intake of sow between the current day and the previous day, while means that the feed intake consumed by the th sow in the current day is smaller than the feed intake consumed in the previous day. However, the same amount of change in the feed intake may reflect some stress level for a sow who typically eats a lot and a more serious stress level for a sow that usually has a lower appetite. For this, we define the relative change in the feed intake by , where is the trimmed average of feed intake of th sow calculated as the average feed intake after removing the lowest and highest of the feed intake measurements taken for the corresponding sow. Here is surrogate for the usual amount of feed intake of the th sow. Trimmed average is used instead of the common average, to remove outliers of very low feed intakes in first few lactating days. For example, consider the situation of two sows: sow that typically consumes 10lb food per day and sow that consumes 5lb food per day. A reduction of 5lb in the feed intake over two consecutive days corresponds to for the th sow and for the th sow. Clearly both sows are stressed (negative value) but the second sow is stressed more, as its absolute relative change is larger; in view of this we refer to the second response as the stress index. Due to the construction of the two types of responses, the data size varies for lactation days 4 (), 11 (), and 18 (); for the first response, , we have sample sizes of 233, 350, and 278, whereas for the sample sizes are 362, 373, and 336 for the respective lactation days.
In this analysis we center the attention on the effect of the ambient temperature and humidity on the 1st quartile of the proxy stress measures and gain more understanding of the food consumption of sows that are most susceptible to heat stress. While the association between the feed intake of lactating sows and the ambient conditions of the farrowing room has been an active research area for some time, accounting for the temperature daily profile has not been considered yet hitherto. Figure 1 displays the temperature and humidity daily profiles recorded at a frequency of 5-minute window intervals for three different days. Preliminary investigation reveals that temperature is negatively correlated with humidity at each time; this phenomenon is caused because the farm uses cool cell panels and fans to control the ambient temperature. Furthermore, it appears that there is a strong pointwise correlation between temperature and humidity. In view of these observations, in our analysis we consider the daily average of humidity. Exploratory analysis of the feed intake behavior of the sows suggests similarities for the sows with parity greater than older sows (ones who are at their third pregnancy or higher); thus we use a parity indicator instead of the actual parity of the sow. The parity indicator is defined as one, if the th sow has parity one and zero otherwise.
For the analysis we smooth daily temperature measurements of each sow using univariate smoother with cubic regression bases and quadratic penalty; REML is used to estimate smoothing parameter. The smoothed temperature curve for sow ’s th repeated measure is denoted by , and the corresponding daily average humidity is denoted by . Both temperature and average humidity are centered before being used in the analysis.
For convenience we denote the response with by removing the superscript. In this application for fixed , corresponds to in Section 2, and correspond to scalar covariates , and and correspond to functional covariates . We first estimate the conditional distribution of given temperature , average humidity , parity , and interaction . Specifically for each of lactation days of interest ( and ) we create a set of equispaced grid of points between the fifth smallest and fifth largest values of ’s and denote the grids with . Then we create artificial binary responses, , and fit the following model for : where is a smooth intercept, quantifies the smooth effect of young sows, describes the effect of the humidity, and and quantify the effect of the temperature at time as well as the interaction between the temperature at time and average humidity. We model using 20 univariate basis functions, and using five univariate basis functions and and using tensor product of two univariate bases functions (total of 25 functions). Throughout the analysis, cubic B-spline bases are used and REML is used for estimating smoothing parameters. The estimated conditional distribution, denoted by , is monotonized by fitting isotonic regression to ; ten smallest and ten largest and the corresponding values of are removed to avoid boundary effects. By abuse of notation, denotes the resulting monotonized estimated distribution. Finally, we obtain estimated first quartiles, i.e., quantiles at level, by inverting , namely, .
To understand the relationship between the lactating sows feed intake and the thermal condition of the farrowing room, we systematically compare and study the predicted quantiles of two responses at combinations of different values of temperature, humidity, and parity. For each of three lactation days () we consider three values of average humidity (first quartile, median, and third quartile) and two levels of parity ( for older sows and for younger sows). Based on the experimenters’ interest, for the functional covariate we consider seven smooth temperature curves given in Figure 2. Each of these curves is obtained by first calculating pointwise quantiles of temperature at five-minute intervals for a specific level and then smoothing it; we considered quantiles levels , and . In short, for each of three lactation days we obtain the first quartile of two responses for 42 different combinations ( humidity values × 2 parity levels × 7 temperature curves) using the proposed method. To avoid extrapolation we ascertain that (i) there are reasonably many observed measurements at each of the combinations and (ii) bottom 25% of the responses are not dominantly from one of the parity group; see distribution of each response by the parity in Figure 3.
The resulting predicted quantiles are shown in Figure 4. Here we focus our discussion on predicted quantile of at quantile level for lactation day 4 ()—the first plot of the second row in Figure 4. The results suggest that the feed intake of older sows (parity ; grey lines) is less affected by high temperatures than younger sows (black lines); this finding is in agreement with Bloemhof et al. [42]. We also observe that the effects of humidity and temperature on feed intake change are strongly intertwined. For illustration, we focus on lactation day 4 () again for younger sows (black lines). For medium humidity (dashed lines) their feed intake stays pretty constant as temperature increases, while for low and high humidity levels (solid and dotted lines, respectively) it changes with an opposite direction. Specifically when temperature increases, the predicted first quartile of increases for low humidity (solid line) whereas it decreases for high humidity (dotted line). Our results imply that high humidity (dotted line) is related to a negative impact of high temperature on feed intake while low humidity (solid line) alleviates it; and this finding is consistent with a previous study [45]. The analysis result suggests to keep low humidity levels in order to maintain healthy feed intake behavior, when ambient temperature is above th percentile; high humidity levels are desirable for cooler ambient temperature.
Interpretation of the other results is similar. While the effects of covariates on feed intake are less apparent toward the end of lactation period, we still observe similar pattern across all three lactation days. For the 11th day (), the 25th quantile of the feed intake is predicted to decrease when the temperature stays below the th percentile, regardless of humidity level and sows age. However, it starts increasing with low humidity while it continues decreasing with high humidity when the temperature rises above the th percentile. Similarly, for the 18th day () when the temperature rises above the th percentile, the predicted first quartile increases with low humidity while it decreases with high humidity. The effect of temperature on feed intake seems less obvious for lactation days 11 and 18 than for day 4; while the effect may be due to lactation day, it may also be a result of other factors, such as more fluctuation and variability in temperature curves on day 4 than on other two days (see Figure 2). Overall we conclude that high humidity and temperature affect the sows feed intake behavior negatively and young sows (parity one) are more sensitive to heat stress than older sows (higher parity), especially at the beginning of lactation period.
6. Discussion
The proposed modeling framework opens up a couple of future research directions. A first research avenue is to develop significance tests of null covariate effect. Testing for the null effect of a covariate on the conditional distribution of the response is equivalent to testing that the corresponding regression coefficient function is equal to zero in the associated function-on-function mean regression model. Such significance tests have been studied when the functional response is continuous [30, 46]; however their study for binary-valued functional responses is an open problem in functional data literature and only recently has been considered in Chen et al. [47]. Another research avenue is to do variable selection in the setting where there are many scalar covariates and functional covariates. Many current applications collect data with increasing number of mixed covariates and selecting the ones that have an effect on the conditional distribution of the response is very important. This problem is an active research area in functional mean regression where the response is normal [48, 49]. The proposed modeling framework has the potential to facilitate studying such problem.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interests.
Acknowledgments
The data used originated from work supported in part by the North Carolina Agricultural Foundation, Raleigh, NC. The authors acknowledge Zhen Han for preparing simulations. Staicu’s research was supported by National Science Foundation DMS 0454942 and DMS 1454942 and National Institutes of Health grants R01 NS085211, R01 MH086633, and 5P01 CA142538-09.
Supplementary Materials
Section S1 provides additional simulation settings and results for the cases of having either a single scalar covariate or a single functional covariate. Additional results for the case of having a scalar covariate and a densely observed functional covariate are also included. Section S2 presents additional data analysis using the proposed method on the bike sharing dataset [50, 51]. (Supplementary Materials)