2 Scalar on composition regression
This case is illustrated using an example from Trinh et al. (2023) studying the impact of climate change on rice yield in Vietnam.
2.1 Rice yield data
The rice_yields dataset contains the rice production
“YIELD” in each province of Vietnam, between 1987 and 2016. The data on
rice yields originate from “International Rice Research Institute”. The
maximum daily temperature data originate from the database “Climate
Prediction Center (CPC)”, created by “National Oceanic and Atmospheric
Administration (NOAA)”. The temperatures have been interpolated over the
63 provinces using areal interpolation methods. For this illustration,
we group the temperatures in three classes LOW ([-6,25.1]),
MEDIUM ([25.1,35.4]) and HIGH ([35.4,45]). The
precipitation data come from the project “Asian
Precipitation-Highly-Resolved Observational Data Integration Towards
Evaluation (APHRODITE)”.
## PROVINCE REGION
## ANGIANG : 30 Central Highlands :180
## BACGIANG : 30 Mekong Delta :420
## BACKAN : 30 Midlands and Northern Mountainous Areas:390
## BACLIEU : 30 Northern and Coastal Central Region :420
## BACNINH : 30 Red River Delta :330
## BARIA-VUNGTAU: 30 Southeastern Area :150
## (Other) :1710
## YEAR YIELD PRECIPITATION
## Min. :1987 Min. :1.120 Min. :0.1693
## 1st Qu.:1994 1st Qu.:3.094 1st Qu.:1.1656
## Median :2002 Median :4.096 Median :1.4040
## Mean :2002 Mean :4.077 Mean :1.3824
## 3rd Qu.:2009 3rd Qu.:5.027 3rd Qu.:1.6411
## Max. :2016 Max. :6.637 Max. :3.7032
##
## TEMPERATURES.LOW TEMPERATURES.MEDIUM TEMPERATURES.HIGH
## Min. :0.0000003 Min. :0.5150684 Min. :0.00000020
## 1st Qu.:0.0027399 1st Qu.:0.6438356 1st Qu.:0.00273983
## Median :0.1945206 Median :0.7657533 Median :0.02185802
## Mean :0.1704263 Mean :0.7957867 Mean :0.03378702
## 3rd Qu.:0.3095891 3rd Qu.:0.9726024 3rd Qu.:0.05205488
## Max. :0.4699453 Max. :0.9999995 Max. :0.25479457
## The following plot is a ternary diagram of maximum temperature colored according to rice yield.
whitered_pal <- colorRampPalette(rev(c("beige","orange","red","darkred")))
top100 <- head(rice_yields, 250)
whitered_col <- whitered_pal(length(top_100$YIELD))## Error in eval(expr, envir, enclos): object 'top_100' not found## Error in eval(expr, envir, enclos): object 'top_100' not found## Error in eval(expr, envir, enclos): object 'whitered_col' not foundlayout(matrix(1:2,ncol=2), width = c(2,1),height = c(1,1))
plot(acomp(top_100$TEMPERATURES),
col = ylcol_top100,
pch = 16)## Error in eval(expr, envir, enclos): object 'top_100' not foundlegend_image <- as.raster(matrix(whitered_pal(20), ncol=1))
plot(c(0,2),c(0,1),type = 'n', axes = F,xlab = '', ylab = '', main = 'Rice yield')
yield_range <- range(top_100$YIELD)## Error in eval(expr, envir, enclos): object 'top_100' not found## Error in eval(expr, envir, enclos): object 'yield_range' not found## Error in eval(expr, envir, enclos): object 'yield_range' not found
2.2 Regression step with the compositions package
This model explains a scalar variable YIELD with a
combination of scalar and compositional variables. Before estimating, we
have to store our data in a data.frame in such a way that scalar
variables are column vectors and compositional variables are matrices.
In our dataset, the TEMPERATURES variable is already in the
matrix format.
## $PROVINCE
## [1] "factor"
##
## $REGION
## [1] "factor"
##
## $YEAR
## [1] "numeric"
##
## $YIELD
## [1] "numeric"
##
## $PRECIPITATION
## [1] "numeric"
##
## $TEMPERATURES
## [1] "matrix" "array"## LOW MEDIUM HIGH
## 0.0000002999999 0.9616434547947 0.0383562452054We continue to fit a regression model with the tools provided by base
R and the compositions package, in particular, the
lm() function and the alr() transformation for
the TEMPERATURES variable.
## [1] "lm"In the next step the function ToSimplex() transforms the
above output into a “lmCoDa” class which serves as the interface between
the lm-output and the interpretation functions.
## [1] "lmCoDa" "lm"The same can be achieved in a single step using the
lmCoDa() function.
fit_X_compo <- lmCoDa(YIELD ~ PRECIPITATION + alr(TEMPERATURES),data = rice_yields)
class(fit_X_compo)## [1] "lmCoDa" "lm"2.2.0.0.1 coef
The estimated parameters can be retrieved by the coef()
function. For a compositional explanatory variable the option
space allows choosing between the simplex version, the clr
version or the coordinate version of the parameters.
## YIELD
## (Intercept) 3.5835270
## PRECIPITATION 0.2528452
## LOW 0.3338069
## MEDIUM 0.3413303
## HIGH 0.3248628## YIELD
## (Intercept) 3.583526963
## PRECIPITATION 0.252845153
## LOW 0.001623843
## MEDIUM 0.023911877
## HIGH -0.025535720## YIELD
## (Intercept) 3.583526963
## PRECIPITATION 0.252845153
## alr(TEMPERATURES)LOW 0.001623843
## alr(TEMPERATURES)MEDIUM 0.0239118772.2.0.0.2 fitted
The fitted() function retrieves the fitted values.
## 1 2 3 4 5 6
## 3.924848 3.736444 3.948922 4.013668 3.871424 3.9522862.3 Interpretation tools
2.3.1 Finite increments interpretation
Finite increments interpretation as in Muller et al. (2018) can be obtained using outputs in coordinate space. For the case of scalar on composition regression, more interpretations are found in Coenders and Pawlovsky-Glahn (2020). The different formulations they propose correspond to different parametrizations of the model and we illustrate below their interpretation based on clr coordinates.
2.3.1.1 Interpretation with ilr coordinates
Using the approach of Muller et al. (2018) we interpret an additive
increment for a specific ilr coordinate that corresponds to a so-called
principal balance. To use their interpretation we first need to
transform the parameters into the appropriate ilr space, which is
defined by the contrast matrix. The ilrBase() function
without further arguments leads to the principal balance of the last
component against all others.
## [,1] [,2]
## 1 -0.7071068 -0.4082483
## 2 0.7071068 -0.4082483
## 3 0.0000000 0.8164966We can use this contrast matrix and the parameters in clr space to obtain the ilr parameters.
clr_TEMP <- coef(fit_X_compo, space = "clr", split = TRUE)[["alr(TEMPERATURES)"]]
ilr_TEMP <- t(V_TEMP) %*% clr_TEMP
ilr_TEMP## YIELD
## [1,] 0.01576002
## [2,] -0.03127474The interpretation is that an additive increment of one unit in the
second ilr coordinate leads to an expected decrease in the rice yield by
0.03 tonnes per hectare. This increment means that the ratio between
HIGH and the geometric mean of LOW and
MEDIUM is approximately multiplied by \(\exp(1/0.82)\sim 3.44,\) where \(O.82\) is the element \((3,2)\) of the contrast matrix.
In terms of temperature distribution, the above increment implies
that the third component HIGH of the temperature
composition grows at the expense of the others (the geometric mean of
the other two) in a way that the ratio of LOW to
MEDIUM remains constant. This corresponds to a change in
the TEMPERATURES simplex in the direction of the vertex
HIGH. Changing for another vertex is simple by a
permutation of the rows of the contrast matrix.
This method only allows to evaluate the impacts of changing a compositional covariate in the direction of a vertex. In Dargel and Thomas-Agnan (2024) it is shown how the impact of more general changes in explanatory compositions can be interpreted.
2.3.1.2 Interpretation with clr coordinates
The clr parameter estimates for the temperatures are given by:
## YIELD
## LOW 0.001623843
## MEDIUM 0.023911877
## HIGH -0.025535720We follow Coenders and Pawlovsky-Glahn (2020) but without changing
the base of the logarithm and the normalization constant. In these
conditions, the first clr parameter for the level LOW of
temperature is interpreted as the expected change of the response when
the ratio between LOW and the geometric mean of
LOW, MEDIUM and HIGH is
multiplied by \(\exp(1)\sim 2.7\). The
resulting change in expected yield is therefore an increase of \(0.0016\) tons per hectare.
2.3.1.3 Variation scenario in the simplex space
To evaluate the impact of a given compositional covariate, we create scenarios of changes in its simplex space described by a linear equation \(x(h)=x(0)\oplus h\odot u,\) where \(x(0)\) is the initial point, \(h\) is the signed-intensity of change and \(u\) is a vector of the simplex defining the direction of change. We then compute the predicted value of \(y(h)\) corresponding to \(x(h)\).
We consider a change in the direction of the vertex “LOW” of the variable “TEMPERATURE”, and we predict new values of the response variable according to the change in percentage points of this vertex.
VariationScenario(
fit_X_compo,
Xvar = "TEMPERATURES",
Xdir = c(0.2, 0.5, 0.3),
inc_size = 2,
n_steps = 5,
add_opposite = TRUE,
obs=43)## YIELD X.LOW X.MEDIUM X.HIGH
## -5 3.819922 0.9989636850642 0.0000018748701 0.0010344400657
## -4 3.855516 0.9963715150059 0.0000314423589 0.0035970426351
## -3 3.891110 0.9870530721012 0.0005237290466 0.0124231988521
## -2 3.926704 0.9498471187296 0.0084740730707 0.0416788081997
## -1 3.962298 0.7674684242080 0.1151254523609 0.1174061234311
## 0 3.997893 0.2465753931507 0.6219176838356 0.1315069230137
## 1 4.033487 0.0220905961112 0.9368346502406 0.0410747536483
## 2 4.069081 0.0013878377613 0.9896156353605 0.0089965268782
## 3 4.104675 0.0000832426310 0.9980354895354 0.0018812678336
## 4 4.140269 0.0000049585605 0.9996043543177 0.0003906871218
## 5 4.175863 0.0000002949985 0.9999186719961 0.0000810330054For a good visual impression, it is better to use a smaller step size and a larger number of steps. Below is the ternary diagram showing the path in the simplex: the initial point is in black and the other ones are colored according to the value of rice yield.
VS2 <- VariationScenario(
fit_X_compo,
Xvar = "TEMPERATURES",
Xdir = c(0.2, 0.5, 0.3),
inc_size = .1,
n_steps = 40,
add_opposite = TRUE,
obs=43)
yield_color <- whitered_pal(length(VS2$YIELD))
yield_color <- yield_color[rank(VS2$YIELD,ties.method = "first")]
plot(acomp(VS2$X), col = yield_color, pch = 16)
plot(acomp(VS2["0", "X"]), add = TRUE, pch = 16)
2.3.2 Infinitesimal increments interpretation
2.3.2.1 Semi-elasticities
The use of semi-elasticities for the interpretation of compositional
covariates is described in Morais and Thomas-Agnan (2021). For a scalar
on composition model the Impacts() method computes the
semi-elasticities for all compositional covariates.
RY_Impacts <- Impacts(fit_X_compo, Xvar = "TEMPERATURES")
RY_Impacts # here, the semi-elasticity is computed for the variable TEMPERATURES ## YIELD
## LOW 0.001623843
## MEDIUM 0.023911877
## HIGH -0.0255357202.3.2.2 Increments approach
Dargel and Thomas-Agnan (2024) use semi-elasticities to interpret
infinitesimal changes due to variation scenarios. The function
VariationTable computes the impact of incremental changes
in covariates on the response variable. The covariate increments,
corresponding to linear changes in the simplex, are indexed by a
direction vector and a signed intensity parameter inc_size
and are specific to a given observation. Two cases are illustrated below
with a direction pointing to a vertex of the simplex and a general
direction. For the vertex direction, the inc_rate is
related to the inc_size parameter and measures the relative
change of the share corresponding to this vertex.
2.3.2.2.1 Infinitesimal increment: direction pointing to a vertex
The direction pointing to a vertex corresponds to the increments
considered by Muller et al. (2018). For example, we can measure the
change caused by increasing the HIGH share of the
TEMPERATURES variable by \(5\%\), assuming that the ratio between the
other two components remains constant, using the function
VariationTable().
VariationTable(
fit_X_compo,
obs = 1, # indicator of the observation (1 is default)
Xvar = "TEMPERATURES", # covariate for which the impact is computed
Xdir = 'HIGH', # vertex in the covariate simplex
inc_rate = 0.05) # signed intensity ## YIELD
## Initial value 3.924848006
## New value 3.919636938
## Semi elasticity -0.031274742
## Variation in % -0.132771206
## Variation in units -0.005211068The interpretation of this result is as follows: a \(5\%\) change of the HIGH share
assuming the ratio between the LOW and the
MEDIUM share remains constant will result in a \(-0.0052\) unit change in
YIELD.
2.3.2.2.2 Infinitesimal increment: general direction
Let us now consider changing the covariate in the direction \((0.2,0.55,0.25)\), with an increment size of \(0.05\). Note that this direction vector will be normalized before being used within the function, and will later be given as output in attributes.
fit_X_compo.varTab <- VariationTable(
fit_X_compo,
obs = 10, # indicator of the observation
Xvar = "TEMPERATURES", # covariate for which the impact is computed
Xdir = c(0.2,0.55,0.25), # general direction in the covariate simplex
inc_size = 0.05) # signed intensity
fit_X_compo.varTab## YIELD
## Initial value 4.010224763
## New value 4.015157689
## Semi elasticity 0.024601746
## Variation in % 0.123008728
## Variation in units 0.004932926## $names
## [1] "YIELD"
##
## $class
## [1] "data.frame"
##
## $row.names
## [1] "Initial value" "New value" "Semi elasticity"
## [4] "Variation in %" "Variation in units"
##
## $`X(0)`
## LOW MEDIUM HIGH
## "0.0000002999999" "0.9506845561646" "0.0493151438356"
## attr(,"class")
## [1] "acomp"
##
## $`X(h)`
## [1] 0.0000002805874 0.9531436237124 0.0468560957002
##
## $Xdir
## [1] 0.1616218 0.6208904 0.2174878
##
## $inc_size
## [1] 0.05
##
## $inc_rates
## [1] -0.064708144 0.002586628 -0.049863955The move between the temperature share vector of \((0.0000002999999, 0.9506845561646,
0.0493151438356)\) to \((0.0000002805874, 0.9531436237124,
0.0468560957002)\) results in an increase of the
YIELD by \(0.004932926\)
units. Note that the Aitchison distance between these two points is
equal to the increment size \(0.05\).
