Conditional analysis for mixed covariates, with application to feed intake of lactating sows

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 involves 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).


Introduction
Many modern applications routinely collect data on study participants comprising scalar responses and both [vector and data streams] information and the main question of interest is to examine how the covariates affect the response. For example in our data application the aim is to study how the daily temperature or humidity behavior in the farrowing rooms (rooms where piglets are born and nursed by the sow until they are weaned) affect the feed intake of lactating sows during their first 21 lactation days, while accounting for other sow-specific information. A popular approach in these cases is to use a semi-parametric framework and assume that the mixed covariates solely affect the mean response; see Cardot et al. (1999); Ramsay and Silverman (2002); James (2002); Vieu (2006, 2009); Goldsmith et al. (2012); McLean et al. (2014) and others.
While it is still important to study the average feed intake, animal scientists are often more concerned with the left tail of the feed intake distribution 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 (1999); ; St-Pierre et al. (2003) among others. In this paper we relax the mean dependence assumption and consider that the covariates affect the entire conditional distribution of the response. Our primary objective is to develop a modeling framework for a comprehensive study of mixed covariates on scalar response; the sow data application includes daily measured temperature and humidity recorded at five minute intervals (functional type) as well as sow-specific information (vector type).
Quantile regression models the effect of scalar/vector covariates beyond the mean response and has attracted great interest (Koenker and Bassett Jr, 1978;Koenker, 2005).
This approach offers a more comprehensive picture of the effects of the covariates on the response distribution. For pre-specified quantile levels, quantile regression models the conditional quantiles of the response as a function of the observed covariates; these approaches have been extended more recently to ensure non-crossing of quantile functions (Bondell et al., 2010). Quantile regression has been also extended to handle functional covariates. Cardot et al. (2005) discussed quantile regression models by employing a smoothing spline modeling-based approach. Kato et al. (2012) 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). Ferraty et al. (2005) and Chen and Müller (2012) considered a different perspective and studied the effect of a functional predictor on the quantiles of the response by first positing a model for the conditional distribution of the response and then inverting it; this approach is appealing as it jointly estimate quantiles for all the desired levels.
More recently, Tang and Cheng (2014), Lu et al. (2014), and Yu et al. (2015) studied quantile regression when the covariates are of mixed types and introduced the partial functional linear quantile regression model framework. The first two papers used fPC basis while the last one proposed to use partial quantile regression (PQR) basis. While such approach is helpful when we are interested in studying the effect of covariate at particular quantile level, it provides an incomplete picture if we are interested in the effect at several quantile levels due to the well-known crossing-issue.
In this paper we propose a comprehensive study of the effect of vector and functional covariates on the distribution of the response. Our approach is inspired from Chen and Müller (2012) (CM, henceforth). Specifically let Q Y |X (τ |X) denote the τ th conditional quantile of Y given a functional covariate X(·), and let F Y |X (y) denote the conditional distribution of Y given X. Using the relationship between Q Y |X (τ |X) and F Y |X (y) that by positing a mean regression model for an auxiliary variable Z(y) := 1(Y ≤ y) and the functional covariate X(·); and 2) estimate Q Y |X (τ |X) by inverting the estimated conditional distribution function. Their estimation approach is restrictive to one functional covariate and a direct extension to accommodate mixed covariates is computationally expensive. We consider a similar idea and propose a modeling framework and estimation technique that easily accommodate various types of covariates in a computationally efficient manner. This development represents the main contribution of our manuscript. Let X 1 be a scalar covariate and X 2 (·) be a functional covariate defined on a closed domain T . We propose the following model for the conditional distribution of Y given X 1 and X 2 (·): where g(·) is a known, monotone link function, β 1 (y) is unknown and smooth function and β 2 (t, y) is unknown and smooth bi-variate function over y and t. The parameters β 1 (y) and β 2 (·, y) quantify the effect of the covariates X 1 and X 2 (·) respectively onto the distribution of the response.
The remainder of the paper is structured as follows. Section 2 discusses the details of the proposed method and Section 3 performs a thorough simulation study evaluating the performance of the proposed method and competitors. We apply the proposed method to analyze the sow data in Section 4. We conclude the paper with a discussion in Section 5.

Modeling Framework
The proposed modeling and estimation method is discussed first for the case of a scalar covariate in Section 2.1; Section 2.2 considers the extension to the case of a functional covariate and then to the case of mixed covariates; Section 2.3 further extends the method to handle sparse and noisy functional covariates. We briefly discuss the monotonization of the estimated conditional distribution in Section 2.4.

Conditional distribution of the response given scalar covariate
First, we focus on the case of a scalar covariate X. Consider the data {(X i , Y i ) : i = 1, . . . , n}, where X i and Y i are independent realizations of real-valued scalar random variables X and Y , respectively. Define Z i (y) = 1(Y i < y) for y ∈ R, where 1(·) is an indicator function; for each y, we view Z i (y) as a binary-valued random variable that is independent and identically distributed as Z(y) = 1(Y < y). It follows that the conditional Here we propose to model the conditional distribution, F Y |X (y), using a generalized function-on-scalar regression model (Goldsmith et al., 2011) between the 'artificial' functional response Z i (y) and the scalar covariate X i .
Specifically, for each y ∈ R, consider where g(·) is a known, monotonic link function, and β 0 (·) and β 1 (·) are unknown, smooth coefficient functions. Here we use the logit function defined as g(x) = log{x/(1 + x)}. It is noteworty to remark that if the slope parameter β 1 (·) is null then the covariate X has no effect on the distribution of the response Y , which is equivalent to X having no effect on any quantile level of Y .
We model β 0 (y) and β 1 (y) by using pre-specified, truncated univariate basis, {B 0,d (·) : where θ 0,d 's and θ 1,d 's are unknown basis coefficients. Then model (2) can be represented as the following generalized additive model where for convenience we use the notation B x,d (y) = XB 1,d (y). The general idea is to set the basis dimensions κ 0 and κ 1 to be sufficiently large to capture the complexity of the coefficient functions and control the smoothness of the estimator through roughness penalties P 0 (θ 0 ) and P 1 (θ 1 ), where θ l , is a vector of all basis coefficients {θ l,d : l = 1, . . . , D l } for l = 1, 2. This approach of using roughness penalties has been widely used; see, for example, Eilers and Marx (1996); Ruppert (2002); Wood (2003Wood ( , 2006a) among many others.
In the following, we detail the estimation algorithm. Let {y j : j = 1, . . . , J} be a set of equi-spaced points in the range of the response variable, Y i 's. For each i and j, we define Z ij = Z i (y j ) = 1(Y i < y j ); it follows that conditional on X i , the Z ij are independently distributed as Bernoulli distribution with mean (µ ij ), where µ ij is such The basis coefficients, θ 0 and θ 1 , are estimated by maximizing the penalized log likelihood criterion, where L is the likelihood function of data {Z ij : j = 1, · · · , J} i , P 0 (θ 0 ) and P 1 (θ 1 ) are penalties, and λ 0 and λ 1 are smoothing parameters. We use quadratic penalties which penalize the size of the curvature of the estimated coefficient functions. Let P 0 (θ 0 ) = θ T 0 D 0 θ 0 and P 1 (θ 1 ) = θ T 1 D 1 θ 1 , where D 0 and D 1 are κ 0 × κ 0 and κ 1 × κ 1 dimensional matrices based on the basis used (see Wood (2006a) for example; the (s, s ) element of D 0 is B 0,s (y)B 0,s (y)dy) and D 1 is defined similarly. The smoothing parameters λ 0 and λ 1 control the trade-off between the goodness of fit and smoothness of the fit. The smoothing parameters are estimated using restricted maximum likelihood (REML).
The criterion (4) can be viewed as the penalized quasi-likelihood (PQL) of the corresponding generalized linear mixed model where D − 0 is the generalized inverse matrix of D 0 and D − 1 is defined similarly. Wood (2006a) discusses an alternative way to deal with the rank-deficient matrices, D 0 and D 1 , in the context of restricted maximum likelihood (REML) estimation. See also Ivanescu et al. (2015) who uses the mixed model representation of a similar regression model to (5), but with a Gaussian functional response.
Let { θ l,d : d = 1, . . . , κ l } for l = 0, 1 be the estimated basis coefficients. It follows that The τ th conditional quantile are estimated as by inverting the estimated distribution This approach relates the τ th level quantile of the response in a nonlinear manner to the covariate. One should note that the estimated distribution function, F Y |X (y), here is not a monotonic function yet. However in practice one can always obtain F Y |X (y) using a monotonization method as described in Section 2.4 first and then invert the resulting estimated distribution to get the estimated conditional quantiles.

Extension to mixed covariates
The modeling approach discussed in Section 2.1 is quite powerful as it can accommodate covariates of various types, including functional covariates.
Assume that a functional covariate X i (·) is observed. For convenience, we assume that X i (·) is observed at fine grid of points and without noise after all. Instead of generalized function-on-scalar model as is used in Section 2.1, now we use generalized function-onfunction regression model E[Z(y)|X] = g −1 β 0 (y) + X(t)β 1 (t, y)dt , where β 1 (·, ·) is an unknown bi-variate coefficient function. In terms of modeling and estimation, the main difference from Section 2.1 is that the coefficient function, β 1 (t, y), is now bivariate and it requires appropriate pre-specified basis function and corresponding penalty term.
The ideas can be further extended to accommodate multiple covariates, scalar or functional, and varied types of effects. For example, assume that we have one vector covariate, one scalar covariate and one functional covariate. We posit a model that considers that the vector covariate has constant effect on the response while both the scalar and the functional covariates have varying effect on the response. Consider the model: for covariate vector X 1 , scalar covariate X 2 , and functional covariate X 3 (·). Here X 1 is assumed to have the y-invariant linear effect, and X 2 and X 3 are assumed to have the y-variant linear effects. It is easy to see that a null effect, say β 3 (·, ·) ≡ 0 is equivalent to the fact that the corresponding covariate, in this case X 3 (·) has no effect on any quantile levels of the response.
Fitting models given in (2) and (6) can be done by extending the ideas of Ivanescu et al. (2015), which provide general modeling and estimation methods for penalized function-on-function regression for the case of Gaussian functional response and implement their method in R (R Core Team, 2013) (namely, the pffr function in refund package (Crainiceanu et al., 2013)). The extension of the model to the non-Gaussian response has recently been studied and implemented by Scheipl et al. (2015). Based on the existing function we implement the proposed method and provide a wrapper function in R.
Furthermore the proposed method can be easily extended to relax the linearity assumption and allow more flexible model structures. For example, instead of a functional linear model such as F Y |X (y) = g −1 β 0 (y) + X 1 β 1 (y) + X 2 (t)β 2 (t, y)dt , we can model the h 1 (·) and h 2 (·, ·, ·) are unknown univariate and trivariate smooth functions, respectively.
We illustrate the nonlinear model in the simulation study for the case when there is a single scalar covariate and the corresponding results are presented in Section S1.1 of the  (Ng and Maechler, 2007), denoted by COBS .
It is important to emphasize that even in the case of a single functional covariate, our methodology differs from Chen and Müller (2012) (CM) in few directions: 1) Our proposed method is based on modeling the unknown smooth coefficient functions using pre-specified basis function expansion and using penalties to control their roughness. In contrast, CM uses data-driven basis, chooses the number of basis functions through the percentage of explained variance (PVE) of the functional predictors. This key difference allows us to accommodate covariates of different types. 2) Our estimation approach is based on a single step penalized criterion while CM uses pointwise estimation based on the residual sum of square criterion and thus requires fitting multiple generalized regressions. This is an important advantage in terms of computational efficiency.

Extension to sparse and noisy functional covariates
In practice the functional covariates are often observed at irregular times across the units and also are possibly corrupted with measurement errors. In such case, one needs to first smooth and de-noise the trajectories before fitting. When the sampling design of the functional covariate is dense, then the common approach is to take each trajectory and smooth it using spline or local polynomial smoothing, as proposed in Ramsay (2006) and Zhang et al. (2007). When the design is sparse, the smoothing is done by pooling all the subjects and using the method proposed in Yao et al. (2005). A method following Yao (Crainiceanu et al., 2013).
In our simulation study we use fpca.sc for this step irrespective of a sampling design (dense or sparse); alternatively, fpca.face (Xiao et al., 2016) in refund or PACE (Yao et al., 2005) in MATLAB (MATLAB, 2010) can also be used.

Monotonization
While a conditional quantile function is nondecreasing, the resulting estimated quantiles may not be. We consider monotonization as proposed by Chernozhukov et al. (2009). Chernozhukov et al. (2009) showed that in this way the monotonized estimator gives the same or better fit than the original estimator. Two approaches are widely used; one is to monotonize the estimated conditional distribution function F Y |X (y), and the other is to monotonize the estimated conditional quantile function Q Y |X (τ ). We choose the former as we already have F Y |X (y) evaluated at dense grid points y j 's and there is no need to obtain the estimated conditional quantile at fine grid points of the quantile level τ ∈ [0, 1]. We use an isotonic regression model (Barlow et al., 1972) for monotonization as it makes no structural assumption and gives an ordered fit; It fits a nonparametric model with an order restriction. The isotonic regression model is fitted through {(y j , F Y |X (y j )) : j = 1, . . . , J} using the isoreg function in R (R Core Team, 2013). This idea was also employed in Kato et al. (2012).

Simulation study
In this section we evaluate numerically the 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 single functional or single scalar covariate are discussed in the Supplementary Materials, Section S1. We adapt the simulation settings of Chen and Müller (2012) for the cases that involve a functional covariate.
The performance is evaluated on a test set of size 100. Two sampling designs are consid- functions with our approach are estimated as described in Section 2 by first creating an artificial binary response Z i (·) and then fitting a penalized function-on-function regression model and using the logit link function; we use the pffr function (Ivanescu et al., 2015;Scheipl et al., 2015) in the refund package (Crainiceanu et al., 2013) in R (R Core Team, 2013) for binomial responses, denoted by Joint QR.
We compare our method with three alternative approaches: (1) a variant of our proposed approach using pointwise fitting, denoted by Pointwise QR, and hence fitting multiple regression models with binomial link function as implemented by the penalized functional regression pfr of the refund package for generalized scalar responses, developed by Goldsmith et al. (2011); (2) a modified version of the CM method, denoted by Mod CM, that we developed to account for additional scalar covariates, and which estimates pointwise using multiple generalized linear models; (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. (2015) and denoted by PQR. Notice that (1) and (2) account for a varying effect of the covariates on the response distribution, but do not ensure that this effect is smooth.
The R function pfr can incorporate both scalar/vector and functional predictors by adopting a mixed effects model framework. The functional covariates are pre-smoothed by fPC analysis and truncation is done based on the number of fPCs determined by a percentage of explained variance equal to 99% for all estimation methods; pre-smoothing the functional covariates before fitting the regression model has been also considered by Goldsmith et al. (2011) and Ivanescu et al. (2015). The performance is evaluated in terms of mean absolute error (MAE) for quantile levels τ = 0.05, 0.1, 0.25, and 0.5.
Numerical results are provided in Tables 1, 2 and 3. Table 1 gives results for the two settings (normal and mixture) when the functional covariate is observed on dense design and the sample size is n = 100; Table 2 shows the corresponding results for n = 1000. We obtained similar findings for the sparse scenario and hence are not reported for brevity.
Consider first the case when conditional distribution of the response is normal. When sample size is large (n = 1000) the proposed method (Joint QR) yields the best MAE for the SNR and the quantile levels considered. Even with low-moderate sample size (n = 100) the Joint QR remains performing the best for extreme quantiles (τ = 0.05, 0.10 and 0.25) and relatively large noises (σ = 4.33 and 6.13). When sample size is smallmoderate the PQR method also performs very well for the small noise (σ = 0.50) and for middle quantiles (τ = 0.05 or τ = 0.10). The Pointwise QR and Mod CM methods perform similarly, where the Pointwise QR tends to do better for low-moderate sample sizes (n = 100) while the Mod CM tends to do better for larger one (n = 1000). All of the four methods are affected by the level of SNR; the higher it is, the better MAE is.
When the conditional distribution of the response follows the mixture of normals, there is no uniformly best method across quantiles levels or SNR levels whe sample size is large (n = 1000). It seems that all four methods have similar performance with some being the best for some situations while others for other situations. Overall the Joint QR method tends to perform better for extreme quantiles (τ = 0.05 or τ = 0.10) while the other three methods tend to predict better the middle quantiles (τ = 0.25 or τ = 0.50). Other findings are relatively similar to the ones for the normal case. 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. The Supplementary Materials, Section S1.1 discusses the former case and compares Joint QR and Pointwise QR with the linear quantile regression and the nonlinear quantile regression (as implemented by the cobs function in the R package COBS (Ng and Maechler, 2007)) 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 methods with nonlinear modeling of the conditional distribution as discussed in Section 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 computational time; see results displayed in Tables S2 and S3. The results show that the proposed Joint QR is superior to CM both in terms of the prediction accuracy and computation efficiency. In our simulation study we also consider the joint fitting of the model by treating the binary response as normal and use pffr (Ivanescu et al., 2015) with Gaussian link, denoted by Joint QR (G).

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 (Rosero et al., 2016). The study comprises of 480 lactating sows of different parities (i.e. the 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 affects 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 Quiniou and Noblet (1999). 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 (Johnston et al., 1999;. Sows with poor feed intake during lactation continue the subsequent reproductive period with negative energy balance (Black et al., 1993), which acts as negative feedback to prevent the onset of a new reproductive cycle. Second, heat stress reduced farrowing rate (the number of sows that deliver a new litter) and reduced the number of piglets born (Bloemhof et al., 2013); The reduction in reproductivity due to seasonality is estimated to cost 300 million dollars per year for the swine industry (St-Pierre et al., 2003). Economic losses are estimated to increase (Nelson et al., 2009) because very high temperatures are likely to occur more frequently due to global warming (Melillo et al., 2014).
Our primary goal is to understand the thermal needs of the lactating sows for proper feeding behavior during the lactation time. Specifically, 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 times during the lactation period: beginning (lactation day 4), middle (day 11) and end (lactation day 18) and consider two types of responses that are meant to assess the feed intake behavior using the current and the previous lactation day. 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 F I ij be the jth measurement of the feed intake observed for the ith sow and denote by LD ij the lactation day when F I ij is measured; here j = 1, . . . , n i . Most sows are observed for every day within the first 21 lactation days and thus have n i = 21. First define the absolute change in the feed intake between two consecutive days as ∆ (1) i(j+1) = 0 means there was no change in feed intake of sow i between the current day and the previous day, while ∆ (1) i(j+1) < 0 means that the feed intake consumed by the ith 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 T A i is the trimmed average of feed intake of ith sow calculated as the average feed intake after removing the lowest 20% and highest 20% of the feed intake measurements taken for the corresponding sow, F I i1 , . . . , F I in i . Here T A i is surrogate for the usual amount of feed intake of the ith 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 i that typically consumes 10lb food per day and sow i that consumes 5lb food per day. A reduction of 5lb in the feed intake over two consecutive days corresponds to ∆ (2) i(j+1) = −10% for the ith sow and ∆ (2) i (j+1) = −20% for the i th sow. Clearly both sows are stressed (negative value) but the second sow is stressed more, as its absolute value is larger; in view of this we refer to the second response as the stress index. Due to the definition of the two types of responses, the data size varies, so for the first response, ∆ (1) i(j+1) , we have sample sizes of 233, 350, and 278 for lactation days 4 (j = 3), 11 (j = 10), and 18 (j = 17), respectively, whereas for ∆ (2) i(j+1) the sample sizes are 362, 373, and 336 respectively. 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 4 displays the temperature and humidity profiles for three 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 strong pointwise correlation between temperature and humidity. In view of these observations we consider the mean summary of the humidity in our analysis.
Exploratory analysis of the feed intake behavior of the sows suggest similarities for the sows with parity greater than one (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 P i is defined as one, if the ith sow has parity one and zero otherwise. As emphasized throughout the paper, the existing literature on quantile regression is not suitable to incorporate covariates of different types, as it is the case here.
For the analysis we smooth daily temperature measurements of each sow using uni-  variate smoother with 15 cubic regression bases and quadratic penalty; REML is used to estimate smoothing parameter. The smoothed temperature curve for sow i's jth repeated measure is denoted by T ij (t), t ∈ [0, 24). In addition the corresponding daily average humidity is denoted by AH ij . Both temperature and average humidity are centered before being used in the analysis.
In the following we detail estimation procedure. Since the procedure is identical for both responses here we remove superscript in notation and use ∆ ij as our response. We first estimate the conditional distribution of ∆ ij given temperature T ij (t), average humidity AH ij , parity P i , and interaction AH ij · T ij (t). Specifically for each of j = 3, 10 and 17 we create a set of 100 equi-spaced grid of points between the fifth smallest and fifth largest values of ∆ ij 's and denote the grids with D = {d : = 1, . . . , 100}. Then we create artificial binary response, {1 (∆ ij < d ) : = 1, . . . , 100}, and fit the following model for where β 0 (·) is a smooth intercept, β 1 (·) quantifies the smooth effect of young sows, β 2 (·) describes the effect of the humidity, and β 3 (·, t) and β 4 (·, t) quantify the effect of the temperature at time t as well as the interaction between the temperature at time t and average humidity. We model β 0 (·) using 20 univariate basis functions, β 1 (·) and β 2 (·) using five univariate basis functions, β 3 (·, ·) and β 4 (·, ·) using tensor product of two univariate bases functions (total of 25 functions). Throughout analysis cubic B-spline bases are used and REML is used for estimating smoothing parameters. The estimated conditional distribution, denoted by F ij (d), is monotonized by fitting isotonic regression to {(d , F ij (d )) : = 10, . . . , 90} ; ten smallest and ten largest d and the corresponding values of F j (d ) are removed to avoid boundary effects. By abuse of notation we again denote the resulting monotonized estimated distribution with F ij (d). Finally we obtain estimated quantiles at τ = 0.25 level by inverting the estimated distribution function To study and interpret the effect of each covariate we predict quantiles at combinations of different values of covariates and investigate the predicted values against the covariates.
For each of three lactation days (j = 3, 10, 17) we consider three values of average humidity (first quartile, median, and third quartile) and the two levels of parity (0 for older sows and 1 for younger sows). For the functional covariate T ij (t) we create seven smooth temperature curves given in Figure 2 by first obtaining pointwise quantiles of temperature and then smoothing them for each of quantiles levels η = 0. at each of the combinations and 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 ∆ i(j+1) at quantile level τ = 0.25 for lactation day 4 (j = 3) -the first plot of the second row in Figure 4. The results suggest that the feed intake of older sows (parity indicator equal to zero) are less affected by high temperatures than that of younger sows; this finding is in agreement with Bloemhof et al. (2013), who also found that younger sows (parity equal to one) are more sensitive to ambient changes than sows with higher parity. We also observe that humidity plays an important role in the effect of temperature on feed intake change. Similar to a previous study (Bergsma and Hermesch, 2012) our results also imply that high humidity is related to a negative impact of high temperature on feed intake while low humidity alleviates it. On the contrary, when temperature is low, high humidity leads to better feed intake than low humidity. For instance on lactation day 4 (j = 3) regardless of the parity group, when temperature increases the predicted first quartile of ∆  order to maintain healthy feed intake behavior, when ambient temperature is above 60th

Response (2) by Parity
percentile; high humidity levels are desirable for cooler ambient temperature.
The results corresponding to other lactation days can also be interpreted similarly.
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 lactation day 11 (j = 10) we observe that when temperature is above 40th percentile the predicted first quartile starts increasing with low humidity while it continues decreasing with high humidity. Similarly for lactation day 18 (j = 17) when temperature is above 60th 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 other two days (see Figure 2). Overall we conclude that high humidity and temperature affect sows' feed intake behavior negatively and young sows (parity one) are more sensitive to heat stress than older sows (higher parity), especially in the beginning of lactation period.

Discussion
In this paper we proposed a novel framework for a comprehensive study of covariates of mixed types on the conditional distribution of the response. Extensive simulation study showed very good prediction performance in terms of estimating quantiles of various levels. Additionally, the modeling approach leads to computationally efficient estimation algorithm. The proposed method is flexible and easy to implement using existent software, pffr .
This 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 i(j+1) and ∆ (2) i(j+1) for different parities, average humidity, and temperature levels. In each of all six panels, black thick lines correspond to the young sows (P i = 1) and grey thin lines correspond to the old sows (P i = 0). Line types indicate different average humidity levels; solid, dashed, and dotted correspond to low, medium, and high average humidity levels (given by the first quartile, median, and the third quartiles of AH ij ), respectively. The seven grids in x-axis of each panel correspond to the 7 temperature curves given in the respective panel of Figure 2.
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 (Shen and Faraway, 2004;Zhang et al., 2007); however their study for binary-valued functional responses is an open problem in functional data literature. One possible alternative is to construct confidence bands for the corresponding coefficient function, say using bootstrap, and examine whether the confidence band includes zero for its entire domain. 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 (Gertheiss et al., 2013;Chen et al., 2015). The proposed modeling framework has the potential to facilitate studying such problem.

Supplementary Material
Section S1 provides additional simulation settings and results for the cases of having either a single scalar covariate or a single functional covariate. Section S2 presents additional data analysis done using the proposed method on the bike sharing dataset (Fanaee-T and Gama, 2013;Lichman, 2013). Lastly the R function implementing the proposed method is available at http://www4.ncsu.edu/~spark13/software/QRFD_Rcode.zip/.