Distributed-lag Linear Structural Equation Models in R: The dlsem Package

In this paper, an extension of linear Markovian structural causal models is introduced, called distributed-lag linear structural equation models (DLSEMs), where each factor of the joint probability distribution is a distributed-lag linear regression with constrained lag shapes. DLSEMs account for temporal delays in the dependence relationships among the variables and allow to assess dynamic causal eﬀects. As such, they represent a suitable methodology to investigate the eﬀect of an external impulse on a multidimensional system through time. In this paper, we present the dlsem package for R implementing inference functionalities for DLSEMs. The use of the package is illustrated through an example on simulated data and a real-world application aiming at assessing the impact of agricultural research expenditure on multiple dimensions in Europe.


Introduction
Structural causal models (SCMs, Pearl 2000, Chapter 5) represent one of the prevalent methodologies for causal inference in contemporary applied sciences.In particular, a Markovian SCM is such that a directed acyclic graph (DAG) encodes causal relationships among the variables, which also implies a factorization of the joint probability distribution according to conditional independence relationships.In a linear parametric formulation (linear Markovian SCM), each factor of the joint probability distribution is a linear regression model, and a causal effect is associated to each edge, directed path or couple of nodes in the DAG to represent average changes in the value of a variable induced by an intervention provoking a unit variation in the value of another variable.A linear Markovian SCM can be extended by letting each factor of the joint probability distribution be a dynamic linear model, in order to account for temporal variations in the dependence relationships among the variables.
When the objective is forecasting, a large variety of dynamic linear models is available.For example, linear regression with polynomial lag shapes (Almon 1965) is a common feature of dynamic predictive models.Mixed-data sampling (MIDAS) regression (Andreou, Ghysels, and Kourtellos 2007;Ghysels, Sinko, and Valkanov 2007), implemented in the R package midasr (Ghysels, Kvedaras, and Zemlys 2016), is a broader class of time series models dealing with data sampled at different frequencies.The R package dlnm (Gasparrini 2011) contains functionalities for linear and non-linear regression models with spline lag-shapes.However, the assessment of causal effects in a dynamic context requires to take into account prior knowledge on lag shapes, and a predictive model may not be adequate without introducing appropriate mathematical constraints.For instance, the effect of an intervention has typically the same sign, or at least is null, at all the possible time lags, but neither the Almon's nor a spline lag shape can be directly applied to this purpose, as they cannot avoid regression coefficients to have different signs.
Constrained lag shapes (Judge, Griffiths, Hill, Lutkepohl, and Lee 1985, Chapters 9-10) overcome such problem, because regression coefficients are non-zero only outside a pre-specified interval of time lags, and within that interval they are allowed to rise from value zero to a maximum before declining again to zero, or to simply decrease from a maximum value to zero.Also, parameter estimation is straightforward, because, given the interval of time lags, the maximum value can be estimated using ordinary least squares.Thus, if prior knowledge on the interval is available, it can be directly exploited, while if the interval is not known a-priori, it may be determined by model selection across all the possible ones.
In this paper, we introduce an extension of linear Markovian SCMs, called distributed-lag linear structural equation models (DLSEMs), where each factor of the joint probability distribution is a distributed-lag linear regression with constrained lag shapes.They were introduced for the first time by Magrini (2018) in the context of lag exposure assessment.DLSEMs represent a suitable methodology to investigate the effect of an external impulse on a multidimensional system through time.
The aim of this paper is to illustrate the dlsem package, which implements inference functionalities for DLSEMs in R. The paper is structured as follows.In Section 2, theory on DLSEMs is presented.In Section 3, instructions for the installation of the dlsem package are provided.In Section 4, the practical use of the dlsem package is illustrated through an example on simulated data.In Section 5, the dlsem package is applied to address a real-world application aiming at assessing the impact of agricultural research expenditure on multiple dimensions in Europe.Section 6 includes conclusive remarks and considerations on future development.

Theory
In this section, theory on distributed-lag linear regression and on structural causal models is provided (Subsections 2.1 and 2.2), then distributed-lag linear structural equation models are presented (Subsection 2.3).

Distributed-lag linear regression
Lagged instances of one or more covariates may be included in the linear regression model to account for temporal delays in their influence on the response: where y t is the value of the response variable at time t and x j,t−l is the value of the j-th covariate at l time lags before t.The set (β j,0 , β j,1 , . . ., β j,L j ) is denoted as the lag shape of the j-th covariate and represents its regression coefficient (in the remainder, simply 'coefficient') at different time lags.
Parameter estimation is inefficient because lagged instances of the same covariate are typically highly correlated.The Almon's polynomial lag shape (Almon 1965) is a well-known solution to this problem, where coefficients for lagged instances of a covariate are forced to follow a polynomial of order Q: Unfortunately, the Almon's polynomial lag shape may show multiple modes and coefficients with different signs, thus entailing problems of interpretation.Constrained lag shapes (Judge et al. 1985, Chapters 9-10) overcome this deficiency.Some examples are the endpointconstrained quadratic lag shape: the quadratic decreasing lag shape: and the gamma lag shape: The endpoint-constrained quadratic lag shape is zero for a lag l ≤ a j − 1 or l ≥ b j + 1, and symmetric with mode equal to θ j at lag (a j + b j )/2.The quadratic decreasing lag shape decreases from value θ j at lag a j to value 0 at lag b j + 1 according to a quadratic function.
The gamma lag shape is positively skewed with mode equal to θ j at lag δ j (δ j −1)log(λ j ) .Value a j is denoted as the gestation lag, value b j as the lead lag, and value b j − a j as the lag width.A static coefficient (no lag shape) is obtained if a j = b j = 0. Since it is not expressed as a function of a j and b j , the gamma lag shape cannot reduce to a static coefficient, but the corresponding values of a j and b j may be computed through numerical approximation from the values of δ j and λ j .
For these three lag shapes it holds: and we refer to the lag sign as the sign of parameter θ j .
A linear regression model with constrained lag shapes is linear in parameters β 0 , θ 1 , . . ., θ J , provided that the values of a 1 , . . ., a J , b 1 , . . ., b J are known.Thus, one may use ordinary least squares to estimate parameters β 0 , θ 1 , . . ., θ J for several models with different values of a 1 , . . ., a J , b 1 , . . ., b J , and then select the one with the minimum value of the Bayesian Information Criterion (BIC, Schwarz 1978) 1 .A heuristic algorithm is shown below.
0. Let J be the number of covariates and T the greatest time lag under consideration.Initialize selected as an empty vector.Set candidates = {X 1 , . . ., X J }.For j = 1, . . ., J, set a j = 0 and b j = 0.
1. Repeat until candidates is not empty: a) initialize fitting as an empty vector.Set cand.temp= candidates; 1 Alternatively, the Akaike Information Criterion (AIC, Akaike 1974) may be used.
b) repeat until cand.temp is not empty: b1) determine j such that the first element in cand.temp is X j .Initialize fit.temp as an empty matrix with T + 1 rows and T + 1 columns; b2) for s 1 = 0, . . ., T and for s 2 = 0, . . ., T : b2.1) set a j = s 1 , b j = s 2 and cov.temp = {selected ∪ X j }; b2.2) use ordinary least squares to estimate parameters β 0 , θ 1 , . . ., θ J for the model with covariates cov.temp, gestation lags {a k : X k ∈ cov.temp} and lead lags {b k : X k ∈ cov.temp}.Compute the BIC for this model and insert it into fit.temp(s 1 , s 2 ); b3) determine k 1 and k 2 such that fit.temp(k 1 , k 2 ) is the best value in fit.temp.
Insert fit.temp(k 1 , k 2 ) into fitting.Set a j = k 1 and b j = k 2 .Remove X j from cand.temp; c) determine X k such that the k-th value in fitting is the minimum one.Insert X k into selected and remove X k from candidates.
Note that neither the response variable nor the covariates must contain unit root in order to obtain unbiased estimates with ordinary least squares (Granger and Newbold 1974).A reasonable procedure is to sequentially apply differentiation to all variables until the Augmented Dickey-Fuller test (Dickey and Fuller 1981) rejects the hypothesis of unit root for all of them.

Structural causal models
Structural causal models (SCMs) were developed by Pearl (2000) in the context of causal inference.They are rooted to path analysis (Wright 1934) and simultaneous equation models (Haavelmo 1943;Koopmans, Rubin, and Leipnik 1950).A SCM consists of a tuple {V , U , Ω V , Ω U , f , P U }, where: -V = {V 1 , . . ., V J } is a set of endogenous variables; - Markovian SCMs (Pearl 2000, Chapter 3) are a special case where f is acyclic and variables in U are each other independent.In a Markovian SCM, the following factorization of the joint probability distribution of variables in V holds: where Π j is the set of variables in V such that, for j > 1, V j is independent of variables in {V 1 , . . ., V j−1 } \ Π j , given variables in Π j .This means that the joint probability distribution of variables in V can be factored according to conditional independence relationships holding among them disregarding variables in U .Pearl (2000, pages 12 and following) shows that these conditional independence relationships are encoded into a directed acyclic graph (DAG) such that Π j is the parent set of V j , ∀ j = 1, . . ., J.For example, in the Markovian SCM associated to the DAG in Figure 1, it holds: and, by way of illustration, V 6 is independent of V 1 and V 2 given V 3 , V 4 and V 5 .
Let do(V i = v i ) denote an intervention setting the value of V i to v i .Then, in a Markovian SCM it holds: where This formula, called truncated factorization (Pearl 2000, Section 3.2), allows to compute the effect of an intervention from the (pre-intervention) distribution in Formula 7, that is to predict the effect of an intervention from non-experimental (observational) data.In a Markovian SCM, the effect of do(V i = v i ) on V j , called causal effect of V i on V j , is given by the following expression (see Pearl 2000, page 70 and following): where Π i is the parent set of V i .
In a linear parametric formulation (linear Markovian SCMs), each factor p(v j | π j ) of the joint probability distribution in Formula 7 is the linear regression model where V j is the response variable and variables in Π j are the covariates.For example, in the linear Markovian SCM associated to the DAG in Figure 1, p(v 4 | v 2 , v 3 ) is the linear regression model where V 4 is the response variable and V 2 and V 3 are the covariates.
In a linear Markovian SCM, the computation of causal effects involves the coefficients of the regression models only, without the need of Formula 10, as shown in the following paragraphs.

Direct causal effects
The coefficient of V i in the regression model of V j , say β j|i , represents the change in the expected value of V j given a unit variation of V i , at constant values of the parents of V j besides V i : Expression 11 is a special case of Expression 10, where the intervention is do(∆V i = 1) and the conditioning set is {Π j \ V i } instead of Π i .Since variables in Π i but not in Π j are independent of V j conditionally to variables in Π j (see Formula 7), we can conclude that β j|i represents the average effect of do(∆V i = 1) on V j : emphasizes that the causal effect in Formula 12 is associated to the edge < V i , V j >.For example, in the linear Markovian SCM associated to the DAG in Figure 1, β 4|3 represents the change in the expected value of V 4 given a unit variation of V 3 , at constant value of V 2 , equating the direct causal effect of V 3 on V 4 .
Indirect causal effects and the overall causal effect Suppose that there exists more than one directed path connecting variable V i to variable V j .In this case, it is straightforward to show that the intervention do(∆V i = 1) influences the expected value of V j independently through each directed path connecting V i to V j , for an overall causal effect equal to the sum of the causal effects associated to each of these paths: where A pathwise causal effect associated to an edge (direct causal effect) can be computed using Formula 12. Instead, a pathwise causal effect associated to a multi-edge directed path, also referred as indirect causal effect, can be computed through the product of the regression coefficients associated to each edge in the path (see, for example, Wright 1934): (14) Note that Formula 2.2 is a generalization of Formula 12.In this view, it is clear that both direct and indirect causal effects belong to the class of pathwise causal effects.For example, in the linear Markovian SCM associated to the DAG in Figure 1, there are three directed paths connecting

Distributed-lag linear structural equation models
Delayed dependence relationships among the variables may be taken into account in a Markovian SCM by specifying each factor of the joint probability distribution in Formula 7 equal to the distributed-lag linear regression in Formula 1.We refer to this Markovian SCM as distributed-lag linear structural equation model (DLSEM).The definition of causal effects at different time lags in a DLSEM is provided in the following paragraphs.

Direct causal effects Let b j (l |
) i be the coefficient of V i at lag l in the regression model of V j .Such coefficient equates the direct causal effect of V i on V j at lag l: l) be the set of all the possible ordered muples of time lags such that their sum is equal to l.If we compute the m direct causal effects associated to each edge in < V d 0 , . . ., V dm > at one of the m-uples in Q ( m l) , say (q 1 , . . ., q m ), and multiply them each other: we obtain one of the possible causal effects of V i on V j through < V d 0 , . . ., V dm > at lag l.Thus, the indirect causal effect of V i on V j through < V d 0 , . . ., V dm > (d 0 = i and d m = j) at lag l is equal to the sum of all the causal effects that can be obtained from Formula 16: Overall causal effects The overall causal effect of V i on V j at lag l, say ∆E (l) (V j | do(∆V i = 1)), is represented by the sum of the pathwise causal effects at lag l associated to each directed path connecting V i to V j .
The cumulative causal effect at a pre-specified time lag, say l, is obtained by summing all the causal effects at each time lag up to l.A pathwise causal lag shape is the set of causal effects associated to a path at different time lags.An overall causal lag shape is the set of the overall causal effects of a variable on another one at different time lags.
The DAG of a DLSEM includes all the possible temporal instances of each variable in V , but it may be represented in a static version for more clarity.For example, only a single temporal instance for each variable is represented, and an edge < V i , V j > exists if and only if there exists at least one time lag where the coefficient of variable V i in the regression model of variable V j is non-zero.
A DLSEM is a special case of dynamic Bayesian network (Murphy 2002), where endogenous variables follow the Gaussian distribution and lag shapes are possibly constrained to predefined functional forms.

Installation
Before installing dlsem, you must have installed R version 2.1.0or higher, which is freely available at http://www.r-project.org/.
To install the dlsem package, type the following in the R command prompt: > install.packages("dlsem")and R will automatically install the package to your system from CRAN.In order to keep your copy of dlsem up to date, use the command: All the results shown in this paper are obtained using R 3.4.3with the dlsem package version 2.3.The official web page of dlsem is: https://cran.r-project.org/web/packages/dlsem/.

Specification of the model code
The first step to build a DLSEM with the dlsem package is the definition of the model code, which includes the formal specification of the regression models.The variables for which a regression model is specified are called endogenous variables.The other variables are referred as exogenous variables.
The model code must be a list of formulas, one for each regression model.In each formula, the response and the covariates must be quantitative variables2 , and operators quec.lag(•),qdec.lag(•) and gamma.lag(•)can be employed to specify, respectively, an endpoint-constrained quadratic, a quadratic decreasing or a gamma lag shape.Operators quec.lag(•) and qdec.lag(•)have three mandatory arguments: • the name of the covariate to which the lag shape is applied, • the gestation lag (a j ), • the lead lag (b j ).
Operator gamma.lag(•)has three mandatory arguments: • the name of the covariate to which the lag shape is applied, • parameter δ j , • parameter λ j .
If none of these two operators is applied to a covariate, it is assumed that its coefficient is equal to 0 for time lags greater than 0 (no lag shape).The group factor and exogenous variables must not appear in the model code (see Subsection 4.3 for the way to include them).
The following code specifies all lag shapes as endpoint-constrained quadratic lag shapes between 0 and 15 time lags: The formula of regression models with no endogenous covariates may be omitted from the model code.For example, the following code (where the formula specifying the regression model for variable Job is omitted) is equivalent to the previous one: > indus.code<-list( + Consum ~quec.lag(Job,0,15),+ Pollution ~quec.lag(Job,0,15)+quec.lag(Consum,0,15) + )

Specification of control options
The second step to build a DLSEM with the dlsem package is the specification of control options, which are distinguished into global (applied to all regression models) and local (model-specific) options.
Global control options must be a named list with one or more of the following components: • adapt: a logical value indicating if adaptation of lag shapes must be performed, that is parameters of lag shapes must be chosen on the basis of fit to data.Default is FALSE, meaning no adaptation; • max.gestation: the maximum gestation lag for all lag shapes.If not provided, it is taken as equal to max.lead (see below); • max.lead: the maximum lead lag for all lag shapes.If not provided, it is computed accordingly to the sample size; • min.width: the minimum lag width for all lag shapes.It cannot be greater than max.lead.If not provided, it is taken as 0; • sign: the lag sign for all lag shapes, that can be either '+' for positive or '-' for negative.If not provided, adaptation will disregard the lag sign; • selection: the criterion to be used for the adaptation of lag shapes, that can be one among 'bic' (the default) and 'aic' to minimise BIC or AIC, respectively.
Local control options must be a named list containing one or more among the following components: • adapt: a named vector of logical values, where each component must have the name of one endogenous variable and indicate if adaptation of lag shapes must be performed for the regression model of that variable; • max.gestation: a named list.Each component of the list must have the name of one endogenous variable and be a named vector.Each component of the named vector must have the name of one covariate in the regression model of the endogenous variable above and include the maximum gestation lag for its lag shape; • max.lead: the same as max.gestation, with the exception that the named vector must include the maximum lead lag; • min.width: the same as max.gestation, with the exception that the named vector must include the minimum lag width; • sign: the same as max.gestation, with the exception that the named vector must include the lag sign (either '+' for positive or '-' for negative).
Local control options have no default values, and global ones are applied in their absence.If some local control options conflict with global ones, only the former are applied.
Suppose that one wants to perform adaptation with the following constraints for all lag shapes: (i) maximum gestation lag of 3 years, (ii) maximum lead lag of 15 years, (iii) minimum lag width of 5 years, (iv) positive lag sign.Control options for these constraints can be expressed in several ways.The most simple solution is to specify only global control options, as the constraints hold for all regression models:

Parameter estimation
Once the model code and control options are specified, parameter estimation can be performed using the command dlsem(•).The main arguments of the command include: • model.code, the first argument of the command, requiring the model code in the format defined in Subsection 4.1; • group, accepting the name of a single group factor (optional).By indicating the group factor, one intercept for each level of the group factor will be estimated in each regression model.
• exogenous, accepting the name of exogenous variables (optional).Exogenous variables can be either qualitative or quantitative, and will be included in each regression model with no lag shape; • log, a logical value indicating whether logarithmic transformation must be applied to all strictly positive quantitative variables (default is FALSE).If log is set to TRUE, the coefficient of each strictly positive quantitative covariate is (approximatively) interpreted as an elasticity, that is as an expected percentage increase in the value of the response variable for 1% increase in the value of the covariate 3 ; • data, requiring an object of class data.framecontaining data; • Messages inform that differentiation was applied one time in order to achieve stationarity of all time series.No mention to imputation was made, meaning that data are complete.The results of command dlsem(•) is an object of class dlsem.Among the components of dlsem objects, we highlight: • estimate: a list of objects of class lm, one for each endogenous variable; 3 The true expected growth rate for the response variable due to 1% increase in the value of a covariate with coefficient κ is equal to 1.01 κ , which corresponds to a percentage increase equal to (1.01 κ − 1) • 100.The approximation (1.01 κ − 1) • 100 ≈ κ here proposed is reasonable for |κ| < 10.
4 Imputation of missing values is performed after eventual logarithmic transformation and differentiation by assuming group-specific means and time-invariant covariance matrix.Qualitative variables cannot contain missing values.Each quantitative variable must have at least 3 observed values if the group factor is not specified, otherwise it must have at least 3 observed values per group.
• model.code: the model code after eventual adaptation; • data.used:data after eventual logarithmic transformation and differentiation.
The summary method for class dlsem returns the summary of the estimation: We see that the number of job positions in industry (Job) significantly influences, on one hand, the amount of private consumption (Consum) from 0 to 5 years and, on the other hand, the amount of greenhouse gas emissions (Pollution) from 1 to 8 years, while the amount of private consumption (Consum) significantly influences the amount of greenhouse gas emissions (Pollution) from 1 to 6 years.This result provides evidence that the influence of industrial development on pollution is both direct and mediated by private consumption.
The plot method for class dlsem displays the DAG of the model in the static representation, where each edge is coloured with respect to the sign of the estimated causal effect (green: positive, red: negative, light gray: not statistically significant): > plot(indus.mod)Note that the DAG includes only the endogenous variables.Argument conf in the plot method controls the confidence level, which is equal to 0.95 by default.Here, a statistically significant and positive causal effect is associated to each edge, thus all of them are shown in green (Figure 3).Colours can be suppressed by setting option style to 1 in the plot method (default is style=2).Instead, by setting option style to 0, all edges are shown in black disregarding statistical significance of causal effects (see Figure 2).

Assessment of causal effects
After parameter estimation is performed by means of command dlsem(•), the command causalEff(•) can be used on the resulting object of class dlsem to compute all the path-

Consum Pollution
Figure 3: The static representation of the DAG for the industrial development problem, where each edge is coloured with respect to the sign of the estimated causal effect.Green: positive causal effect.Red: negative causal effect.Light grey: not statistically significant causal effect (no such edges here).
wise causal lag shapes and the overall one connecting two variables.The main arguments of command causalEff(•) include: • x, the first argument, requiring an object of class dlsem; • from, requiring the name of one or more starting variables, that is the variables generating the causal effect; • to, requiring the name of the ending variable, that is the variable receiving the causal effect; • lag, accepting specific time lags at which the causal effect must be computed.If no values are provided, all the relevant time lags are considered.
• cumul, a logical value indicating whether the cumulative causal effect must be returned.Default is FALSE.
Only exogenous variables can be indicated as starting or ending variables.Note that, due to the properties of the multiple linear regression model, causal effects are net of the influence of the group factor and exogenous variables.
Since the logarithmic transformation was applied to all quantitative variables, the resulting causal effects are interpreted as elasticities, that is, for a 1% of job positions more, greenhouse gas emissions are expected to grow by 0.61% after 5 years and by 1.11% after 10 years.The influence ends after 11 years, as the cumulative causal effects at lags 11 and 12 are equal.

Pollution ~ Job
Lag Coefficient q q q q q q q q q q q q q 0 1 2 > lagPlot(indus.mod,from="Job",to="Pollution") The resulting graphics are shown in Figure 4.Note that a multi-edge pathwise causal lag shape is a mixture of different lag shapes, thus it may show an irregular aspect, like it is the case of the overall causal lag shape displayed in the lower panel of Figure ??.

Model comparison
We now fit two alternative models for the industrial development problem, such that all lag shapes are quadratic decreasing and gamma lag shapes, respectively: We see that the three models provide different results.The BIC method for class dlsem can be used to compare them according to the BIC: > lapply(list(QUEC=indus.mod,QDEC=indus.mod_2,GAMMA=indus.mod_3),BIC) $QUEC Job Consum Pollution (overall) -1733.060 -1553.811 -1349.505 -4636.377$QDEC Job Consum Pollution (overall) For further details, see the documentation of the agres dataset by typing ?agres.
In order to highlight the opportunity to include qualitative exogenous variables, we consider an additional variable in the context level, that is a dummy indicating whether the Decoupling policy implemented in 2005 is in vigour or not.Such variable can be defined with the following code: > agres$POLICY <-factor(1*(agres$YEAR>=2005)) > levels(agres$POLICY) <-c("no","yes") It is important to indicate qualitative variables with numerical labels as factors, otherwise they will be interpreted as quantitative variables.With the following code we define the four causal levels and request the summary of the considered data: > context.var<-c("GDP","EMPL_AGR","UAA","PATENT_OTHER","POLICY") > investment.var<-c("GBAORD_AGR","BERD_AGR") > research.var<-c("RD_EDU_AGR","PATENT_AGR") > impact.var<-c("ENTR_INCOME_AGR","PPI_AGR") > all.var<-c(context.var,investment.var,research.var,impact.varFrom these summaries it seems that propensity to research investments is heterogeneous across countries: those with higher gross domestic product and utilized agricultural area tend to expend more for agricultural research.The summary of parameter estimation can be obtained using the summary method.Here, we display the summary for the endogenous part only by accessing to the component endogenous: The static representation of the DAG with coloured edges is shown in Figure 6: > plot(agres.mod) Results show that business enterprise research expenditure (BERD_AGR) influences producer price (PPI_AGR) both directly from 3 to 20 years, and indirectly through the number of researchers with tertiary education (RD_EDU_AGR) from 0 to 5 years.Producer price (PPI_AGR) is also influenced by the number of patent applications (PATENT_AGR) from 1 to 6 years, independently of business enterprise research expenditure (BERD_AGR), which also influences entrepreneurial income (ENTR_INCOME_AGR) from 0 to 12 years.

RD_EDU_AGR PATENT_AGR ENTR_INCOME_ PPI_AGR
Figure 6: The static representation of the DAG for the agricultural research problem (endpoint-constrained lag shapes) where each edge is coloured with respect to the sign of the estimated causal effect.Green, red and light grey indicate positive, negative and not statistically significant causal effects, respectively.

Lag
Coefficient q q q q q q q q q q q q q q 0 1 2

Lag
Coefficient q q q q q q q q q q q q q q q q q q q q q q q 0 1 2 3 4 5 6 7 8 The static representation of the DAG with coloured edges is shown in Figure 8: According to the BIC, the model with gamma lag shapes should be preferred:

Conclusions and future development
Distributed-lag linear structural equation models (DLSEMs) allow to study the dynamic causal path generated by an external impulse into a multi-dimensional system, thanks to a variety of lag shapes characterized by flexibility and simplicity of estimation.In the real-world application illustrated in this paper, DLSEMs are applied to an impact assessment problem in the field of agricultural economics, but they appear as a suitable solution in several other scientific domains where the effect of an external impulse is studied through time.For instance, epidemiologic problems may be addressed, like the analysis of the dynamic effect of pollutants on human health (Martins, Pereira, Lin, Santos, Prioli, do Carmo Luiz, Saldiva, and Ferreira Braga 2006;Steinvil, Fireman, Kordova-Biezuner, Cohen, Shapira, Berliner, and Rogowski 2009).Also, viticultural-oenological domains can be dynamically characterized by focusing, for example, on grape maturation (Magrini, Di Blasi, and Stefanini 2017) and vinification (Stefanini and Pantani 2013;Magrini, Pantani, Bartolini, and Stefanini 2016).Lag shapes included in the package may represent a large number of real-world lag structures: unimodal symmetric (with the endpoint-constrained quadratic lag shape), unimodal asymmetric (with the gamma lag shape) and skewned ones (with the quadratic decreasing lag shape).Nevertheless, additional lag shapes with further specific features may be added in future.
Parameter estimation in DLSEM cannot be performed in a single step unless gestation and lead lags are all known.Since the number of possible models rises exponentially as the number of covariates and time lags increases, complete search is infeasible for most real-world applications, thus a heuristic search was implemented.Further development of the package may be directed towards the improvement of the search strategy.
The main limitation of DLSEMs relies in the sample size to achieve efficient estimates.Likewise all lag models proposed in the literature, DLSEMs require to drop away a number of statistical units in order to estimate as many regression coefficients at different time lags.Thus, long time series are typically required to efficiently estimate wide lag structures.An extensive simulation study to assess the potential of the method as a function of the sample size is a natural continuation of the research on this topic.
DLSEMs are a special case of dynamic Bayesian networks (Murphy 2002), where endogenous variables follow the Gaussian distribution and lag shapes are possibly constrained to predefined functional forms.An important functionality of DLSEMs implemented in the dlsem package is the opportunity to condition on exogenous variables of any type, either categorical or continuous.
The current implementation of the package deals with grouped data through fixed effects estimation.Feature releases may include random effects estimation to enhance inference whenever the considered groups are a subset of the possible ones, or covariates with values constant within groups (second-level covariates) are involved.

Figure 1 :
Figure 1: An example of directed acyclic graph.

Figure 4 :
Figure 4: The pathwise causal lag shapes (upper panels) and the overall one (lower panel) connecting Job and Pollution.95% asymptotic confidence intervals are shown in grey.

Figure 5 :
Figure 5: The static representation of the DAG for the agricultural research problem.

Figure 8 :
Figure 8: The static representation of the DAG for the agricultural research problem (gamma lag shapes) where each edge is coloured with respect to the sign of the estimated causal effect.Green, red and light grey indicate positive, negative and not statistically significant causal effects, respectively.
The static representation of the DAG for the industrial development problem.'Job': number of job positions in industry.'Consum': private consumption index.'Pollution': amount of greenhouse gas emissions.
Figure 7: The pathwise causal lag shapes (upper panels) and the overall one (lower panel) connecting BERD_AGR to PPI_AGR.95% asymptotic confidence intervals are shown in grey.