\documentclass[a4paper,DIV12]{scrartcl} \usepackage{amsmath} \usepackage{amsfonts} \usepackage{color} \usepackage[bookmarks=TRUE, colorlinks, citecolor=black, filecolor=black, linkcolor=black, urlcolor=black, ]{hyperref} \newcommand{\code}[1]{\texttt{#1}} \newcommand{\se}{\mathrm{se}\,} \newcommand{\var}{\mathrm{Var}\,} \setlength{\emergencystretch}{3em} %\VignetteEngine{knitr::knitr} %\VignetteIndexEntry{Extracting and Unifying Semi-Elasticities and Effect Sizes from Studies with Binary and Categorical Dependent Variables} \title{ Extracting and Unifying Semi-Elasticities\\ and Effect Sizes from Studies\\ with Binary and Categorical Dependent Variables} \author{Geraldine Henningsen \& Arne Henningsen} \begin{document} %\SweaveOpts{concordance=TRUE} \maketitle This is an addendum to our vignette ``Extracting and Unifying Semi-Elasticities and Effect Sizes from Studies with Binary Dependent Variables.'' %In this appendix we will give detailed account of the way we derived %the semi-elasticities from all articles used in the meta-study for the %variables income, age, education, household size, and home-ownership. %Furthermore, the appendix reports the details of all studies used in table X %and in the meta-study, as well, as the R-code used to generate the %semi-elasticities from each study used in the meta-analysis. %The variables of interest, income, age, education, household size, and %home-ownership are standardised as follows: %\begin{itemize} %\item[]\textit{income} - continuous variable %\item[]\textit{age} - categorical variable with three categories: young (ref.) with % mean 27 years, mid-age with mean 45 years, old-age with mean 70 years %\item[]\textit{education} - binary variable with low education (no college/university) % as reference group and high education (college/university) as effect group %\item[]\textit{household size} - continuous variable %\item[]\textit{home-ownership} - binary variable with renter as reference group and ownership % as effect group %\end{itemize} %============================================================================================== \section{Linear probability models and articles reporting marginal effects} \subsection{Semi-elasticities from continuous explanatory variables (linear and quadratic)} \label{sec:LPMcont} Linear probability model with continuous explanatory variables:% \footnote{The following explanations are focused on linear probability models but apply equally to articles with probit or logit regressions that report the marginal effects.} \begin{equation} y = \beta_0 + \sum_{k=1}^K \beta_k x_k + u, \end{equation} where $y \in \{ 0, 1 \}$ is a binary dependent variable, $x = ( x_1, \ldots , x_K )^{\top}$ is a vector of $K$ continuous explanatory variables, $\beta = ( \beta_0 , \ldots , \beta_K )^{\top}$ is a vector of $K + 1$ unknown coefficients, and $u$~is a random error term. The semi-elasticity of the $k$th (continuous) explanatory variable is: \begin{equation} \epsilon_k = \frac{\partial E[y]}{\partial x_k} \; x_k = \beta_k \cdot x_k . \label{eq:linEla} \end{equation} This semi-elasticity can be interpreted as: if the $k$th explanatory variable inceases by one percent, the probability that $y=1$ increases by $\epsilon_k$ percentage points. This semi-elasticity depends on the value of the dependent variable~$x_k$. One often calculates the semi-elasticity at the sample mean, i.e., $x_k = \bar{x}_k$. If the model specification additionally includes a quadratic term of the explanatory variable~$k$, e.g., $x_{k+1} = x_k^2$, the semi-elasticity of this explanatory variable is: \begin{equation} \epsilon_k = \frac{\partial E[y]}{\partial x_k} \; x_k = \left( \beta_k + 2 \; \beta_{k+1} \; x_k \right) x_k = \beta_k \; x_k + 2 \; \beta_{k+1} \; x_k^2 . \label{eq:linElaQuad} \end{equation} % \subsubsection{ Standard errors for semi-elasticities for continuous variables} An approximate standard error of the semi-elasticity defined in equation~(\ref{eq:linElaQuad}) can be obtained by using the Delta-method: \begin{align} \se( \epsilon_k ) & = \sqrt{ \frac{ \partial \epsilon_k }{ \partial ( \beta_k , \beta_{k+1} ) } \var( \beta_k, \beta_{k+1 } ) \frac{ \partial \epsilon_k }{ \partial ( \beta_k , \beta_{k+1} )^{\top} } } \\ & = \sqrt{ \left( x_k , 2 \; x_k^2 \right) \; \var( \beta_k, \beta_{k+1} ) \; \left( x_k , 2 \; x_k^2 \right)^{\top} }, \label{eq:linElaQuadSE} \end{align} where $\se( \epsilon_k )$ indicates the (approximate) standard error of the semi-elasticity $\epsilon_k$ and $\var( \beta_k, \beta_{k+1} )$ indicates the variance-covariance matrix of $\beta_k$ and $\beta_{k+1}$. As scientific publications usually do not report covariances between estimated coefficients, but only (at best) standard errors (or t-values, which can be used to calculate the standard errors), the covariance between the estimates of $\beta_k$ and $\beta_{k+1}$ is usually unknown. One can approximate equation~(\ref{eq:linElaSE}) by assuming that the covariance between the estimates of $\beta_k$ and $\beta_{k+1}$, i.e., the off-diagonal element(s) of $\var( \beta_k, \beta_{k+1} )$, is zero: \begin{equation} \var( \beta_k, \beta_{k+1} ) \approx \left[ \begin{array}{cc} \se( \beta_k )^2 & 0 \\ 0 & \se( \beta_{k+1} )^2 \end{array} \right], \label{eq:linElaQuadVar} \end{equation} where $\se( \beta_k )$ and $\se( \beta_{k+1} )$ are the standard errors of the estimates of $\beta_k$ and $\beta_{k+1}$, respectively. If there is no quadratic term of an explanatory variable~$x_k$ and, thus, the semi-elasticity is calculated according to equation~(\ref{eq:linEla}), equation~(\ref{eq:linElaQuadSE}) simplifies to: \begin{equation} \se( \epsilon_k ) = \se( \beta_k ) \cdot x_k \label{eq:linElaSE} \end{equation} The following function calculates the semi-elasticity~$\epsilon_k$ according to equation~(\ref{eq:linElaQuad}) and its standard error according to equations~(\ref{eq:linElaQuadSE}) and~(\ref{eq:linElaQuadVar}): <<>>= linEla <- function( xCoef, xVal, xCoefSE = rep( NA, length( xCoef ) ) ){ if( ! length( xCoef ) %in% c( 1, 2 ) ) { stop( "argument 'xCoef' must be a scalar or vector with 2 elements" ) } if( length( xCoef ) != length( xCoefSE ) ) { stop( "arguments 'xCoef' and 'xCoefSE' must have the same length" ) } if( length( xVal ) != 1 || !is.numeric( xVal ) ) { stop( "argument 'xVal' must be a single numeric value" ) } if( length( xCoef ) == 1 ) { xCoef <- c( xCoef, 0 ) xCoefSE <- c( xCoefSE, 0 ) } semEla <- ( xCoef[1] + 2 * xCoef[2] * xVal ) * xVal derivCoef <- c( xVal, ifelse( xCoef[2] == 0, 0, 2 * xVal^2 ) ) vcovCoef <- diag( xCoefSE^2 ) semElaSE <- drop( sqrt( t( derivCoef ) %*% vcovCoef %*% derivCoef ) ) result <- c( semEla = semEla, stdEr = semElaSE ) return( result ) } @ with argument \code{xCoef} $= ( \beta_k , \beta_{k+1} )^{\top}$, argument \code{xVal} $= x_k$, and an optional argument \code{xCoefSE} $= ( \se( \beta_k ), \se( \beta_{k+1} ) )^{\top}$. In case of model equations that are linear in $x_k$, the second elements of arguments \code{xCoef} and \code{xCoefSE} can be set to zero or can be omitted, i.e., \code{xCoef} $= \beta_k$ and \code{xCoefSE} $= \se( \beta_k )$, so that equation~(\ref{eq:linElaQuad}) simplifies to equation~(\ref{eq:linEla}) and equation~(\ref{eq:linElaQuadSE}) simplifies to equation~(\ref{eq:linElaSE}). << >>= # Example # model equation that is linear in x_k ela1a <- linEla( 0.05, 23.4 ) ela1a ela1b <- linEla( 0.05, 23.4, 0.001 ) ela1b all.equal( ela1b, linEla( c( 0.05, 0 ), 23.4, c( 0.001, 0 ) ) ) # Example # model equation that is quadratic in x_k ela2a <- linEla( c( 0.05, -0.00002 ), 23.4 ) ela2a ela2b <- linEla( c( 0.05, -0.00002 ), 23.4, c( 0.001, 0.00002 ) ) ela2b @ %======================================================================================= \subsection{Semi-elasticities from interval-coded explanatory variables} \label{sec:linInterval} Linear probability model, where the $k$th explanatory variable is interval-coded: \begin{align} y & = \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m \in \{ 1, \ldots , M \} \setminus m^* } \delta_m D_m + u,\\ D_m & = \begin{cases} 1 & \text{ if } b_{m-1} < x_k \leq b_m \\ 0 & \text{ otherwise} \end{cases} \quad \forall \; m = 1, \ldots , M, \end{align} where $y \in \{ 0, 1 \}$ is a binary dependent variable, $x = ( x_1, \ldots , x_K )^{\top}$ is a vector of $K$ continuous explanatory variables, whereas the actual values of one of these variables, $x_k$, are unobserved, $D = ( D_1 , \ldots , D_M )^{\top}$ is a vector of $M$~dummy variables that indicates in which intervals the values of variable $x_k$ fall, $b = ( b_0 , \ldots , b_M )^{\top}$ is a vector of the $M + 1$ boundaries of the $M$~intervals of variable~$x_k$, $m^* \in \{ 1 , \ldots , M \}$ is an arbitrary chosen interval that is used as `base' interval in the regression, $\beta = ( \beta_0 , \ldots , \beta_K )^{\top}$ and $\delta = ( \delta_1, \ldots, \delta_{m^* - 1}, \delta_{m^* + 1} , \ldots, \delta_M )^{\top}$ are vectors of $K + 1$ and $M - 1$, unknown coefficients, respectively, and $u$~is a random error term. For convenience of further calculations, we define the (non-estimated) coefficient for the `base' interval to be zero, i.e., $\delta_{m^*} = 0$. The semi-elasticity of explanatory variable~$x_k$: \begin{equation} \epsilon_k \equiv \frac{ \partial E[y] }{ \partial x_k } \; x_k \end{equation} can be approximated `around' each `inner' interval boundary $b_1 , \ldots , b_{M-1}$ by: \begin{align} \epsilon_{km} \; \approx & \frac{ E[ y | b_m < x_k \leq b_{m+1} ] - E[ y | b_{m-1} < x_k \leq b_m ] }{ E[ x_k | b_m < x_k \leq b_{m+1} ] - E[ x_k | b_{m-1} < x_k \leq b_m ] } \; b_m \label{eq:linElaIntBound}\\[1ex] = & \frac{ \delta_{m+1} - \delta_m }{ E[ x_k | b_m < x_k \leq b_{m+1} ] - E[ x_k | b_{m-1} < x_k \leq b_m ] } \; b_m \label{eq:linElaIntM} \\[1ex] & \quad \forall \; m = 1, \ldots , M-1 . \nonumber \end{align} It indicates the approximate increase in the probability of $y=1$ (in percentage points) that is caused by an increase in the explanatory variable~$x_k$ by one percent around the interval border~$b_m$. Assuming: \begin{equation} E[ x_k | b_{m-1} < x_k \leq b_m ] \approx \frac{1}{2} \left( b_{m-1} + b_m \right), \end{equation} we get: \begin{align} \epsilon_{km} \; \approx & \frac{ \delta_{m+1} - \delta_m }{ \frac{1}{2} \left( b_m + b_{m+1} \right) - \frac{1}{2} \left( b_{m-1} + b_m \right) } \; b_m \\[1ex] = & 2 \; \frac{ \delta_{m+1} - \delta_m }{ b_{m+1} - b_{m-1} } \; b_m \label{eq:linElaInt} \\[1ex] & \quad \forall \; m = 1, \ldots , M-1 \nonumber \end{align} If $b_M$ is infinity, we can assume: \begin{align} E[ x_k | b_{M-1} < x_k \leq b_M ] & \approx b_{M-1} + \left( b_{M-1} - b_{M-2} \right) \\ & = 2 \; b_{M-1} - b_{M-2} \\ & = \frac{1}{2} \left( 4 \; b_{M-1} - 2 \; b_{M-2} \right) \\ & = \frac{1}{2} \left( b_{M-1} + 3 \; b_{M-1} - 2 \; b_{M-2} \right), \end{align} which is equivalent to setting: \begin{equation} b_M = 3 \; b_{M-1} - 2 \; b_{M-2}. \label{eq:linElaIntBM} \end{equation} In order to aggregate the semi-elasticities $\epsilon_{k1}, \ldots , \epsilon_{k,M-1}$ that correspond to the `inner' interval boundaries $b_1 , \ldots , b_{M-1}$ to one overall semi-elasticity, we can calculate a weighted mean of the semi-elasticities $\epsilon_{k1}, \ldots , \epsilon_{k,M-1}$. We suggest to use the approximate proportions of observations that have a value of variable~$x_k$ that are `around' each boundary $b_1 , \ldots , b_{M-1}$ as weights: \begin{equation} w_m = \begin{cases} s_1 + \frac{1}{2} s_2 & \text{ if } m = 1 \\ \frac{1}{2} \left( s_m + s_{m+1} \right) & \text{ if } 2 \leq m \leq M-2 \\ \frac{1}{2} s_{M-1} + s_M & \text{ if } m = M-1 \end{cases}, \label{eq:linElaIntWeights} \end{equation} where $s_m$ is the proportion of observations that are in the $m$th interval, i.e., $b_{m-1} < x_k \leq b_m$. The weights sum up to one, i.e., $\sum_{m=1}^{M-1} w_m = 1$. Thus, we can calculate the approximate average semi-elasticity by: \begin{align} \epsilon_k \approx \; & \sum_{m=1}^{M-1} w_m \; \epsilon_{km} \label{eq:linElaIntAvg}\\ = & \left( 2 s_1 + s_2 \right) \frac{ \delta_2 - \delta_1 }{ b_2 - b_0 } \; b_1 \label{eq:linElaIntAvgFull}\\ & + \sum_{m=2}^{M-2} \left( s_m + s_{m+1} \right) \frac{ \delta_{m+1} - \delta_m }{ b_{m+1} - b_{m-1} } \; b_m \nonumber \\ & + \left( s_{M-1} + 2 s_M \right) \frac{ \delta_M - \delta_{M-1} }{ b_M - b_{M-2} } \; b_{M-1}. \nonumber \end{align} As argued before, if $b_M$ is infinity, it can be set to the right-hand side of equation~(\ref{eq:linElaIntBM}) or another value that seems to be appropriate for the analysed data set. An approximate standard error of the semi-elasticity defined in equation~(\ref{eq:linElaIntAvgFull}) can be obtained by using the Delta-method: \begin{align} \se( \epsilon_k ) & = \sqrt{ \left( \frac{ \partial \epsilon_k }{ \partial \delta } \right)^{\top} \var( \delta ) \; \frac{ \partial \epsilon_k }{ \partial \delta } }, \label{eq:linElaIntSE} \end{align} where $\se( \epsilon_k )$ indicates the (approximate) standard error of $\epsilon_k$ and $\var( \delta )$ indicates the variance-covariance matrix of the estimates of~$\delta$. The $n$th element of the vector of partial derivatives $\partial \epsilon_k / \partial \delta$ can be obtained by: \begin{align} \frac{ \partial \epsilon_k }{ \partial \delta_n } & = \sum_{m=1}^{M-1} w_m \; \frac{\partial \epsilon_{km}}{\partial \delta_n}\\ & = \begin{cases} w_1 \; \dfrac{\partial \epsilon_{k1}}{\partial \delta_1} & \text{ if } n = 1 \\[2ex] w_{n-1} \; \dfrac{\partial \epsilon_{k,n-1}}{\partial \delta_n} + w_n \dfrac{\partial \epsilon_{kn}}{\partial \delta_n} & \text{ if } 2 \leq n \leq M-1 \\[2ex] w_{M-1} \; \dfrac{\partial \epsilon_{k,M-1}}{\partial \delta_M} & \text{ if } n = M \end{cases} \\ & = \begin{cases} - 2 \; w_1 \; \dfrac{ b_1 }{ b_2 - b_0 } & \text{ if } n = 1 \\[2ex] 2 \; w_{n-1} \; \dfrac{ b_{n-1} }{ b_n - b_{n-2} } - 2 \; w_n \; \dfrac{ b_n }{ b_{n+1} - b_{n-1} } & \text{ if } 2 \leq n \leq M-1 \\[2ex] 2 \; w_{M-1} \; \dfrac{ b_{M-1} }{ b_M - b_{M-2} } & \text{ if } n = M \end{cases} \label{eq:secases} \end{align} The following code defines a helper function that calculates the weights~$w_1, \ldots, w_{M-1}$ according to equation~(\ref{eq:linElaIntWeights}): <<>>= elaIntWeights <- function( xShares ) { nInt <- length( xShares ) weights <- rep( NA, nInt - 1 ) for( m in 1:(nInt-1) ){ weights[m] <- ifelse( m == 1, 1, 0.5 ) * xShares[m] + ifelse( m+1 == nInt, 1, 0.5 ) * xShares[m+1] } if( abs( sum( weights ) - 1 ) > 1e-5 ) { stop( "internal error: weights do not sum up to one" ) } return( weights ) } @ with argument \code{xShares} $= s = ( s_1 , \ldots , s_M )^{\top}$ the vector of proportion of observations in each interval. The following code defines a helper function that checks the boundaries~$b_0, \ldots, b_M$ and replaces $b_M$ by the right-hand side of equation~(\ref{eq:linElaIntBM}) if $b_M$ is infinite: <>= elaIntBounds <- function( xBound, nInt, argName = "xBound" ) { if( length( xBound ) != nInt + 1 ) { stop( "argument '", argName, "' must be a vector with ", nInt + 1, " elements" ) } if( any( xBound != sort( xBound ) ) ) { stop( "the elements of the vector specified by argument '", argName, "' must be in increasing order" ) } if( max( table( xBound ) ) > 1 ) { stop( "the vector specified by argument '", argName, "' may not contain two (or more) elements with the same value" ) } if( is.infinite( xBound[ nInt + 1 ] & nInt > 1 ) ) { xBound[ nInt + 1 ] <- 3 * xBound[ nInt ] - 2 * xBound[ nInt - 1 ] } return( xBound ) } @ with argument \code{xBound} $= b = ( b_0 , \ldots , b_M )^{\top}$ the vector of boundaries, argument \code{nInt} $= M$ an integer that indicates the number of intervals, and optional argument \code{argName} a character string that can be used to modify the error message (if necessary). Using the helper functions \code{elaIntWeights} (equation~\ref{eq:linElaIntWeights}) and \code{elaIntBounds} (equation~\ref{eq:linElaIntBM}), the following code defines a function that calculates the semi-elasticity~$\epsilon_k$ according to equations~(\ref{eq:linElaInt}) and~(\ref{eq:linElaIntAvg}) and the respective standard error according to equations~(\ref{eq:linElaIntSE}) and (\ref{eq:secases}): << >>= linElaInt <- function( xCoef, xShares, xBound, xCoefSE = rep( NA, length( xCoef ) ) ){ nInt <- length( xCoef ) if( nInt < 2 || !is.vector( xCoef ) ) { stop( "argument 'xCoef' must be a vector with at least two elements" ) } if( length( xCoefSE ) != nInt ) { stop( "arguments 'xCoef' and 'xCoefSE' must be vectors of the same length" ) } if( length( xShares ) != nInt ) { stop( "arguments 'xCoef' and 'xShares' must be vectors of the same length" ) } if( any( xShares < 0 ) ) { stop( "all shares in argument 'xShares' must be non-negative" ) } if( abs( sum( xShares ) - 1 ) > 0.015 ) { stop( "the shares in argument 'xShares' must sum to one" ) } # check 'xBound' and replace infinite values xBound <- elaIntBounds( xBound, nInt ) # weights weights <- elaIntWeights( xShares ) # semi-elasticities 'around' each inner boundary and their weights semElas <- rep( NA, nInt - 1 ) for( m in 1:(nInt-1) ){ semElas[m] <- 2 * ( xCoef[ m+1 ] - xCoef[ m ] ) * xBound[ m+1 ] / ( xBound[m+2] - xBound[m] ) } # (average) semi-elasticity semElaAvg <- sum( semElas * weights ) # derivatives of the (average) semi-elasticity wrt the coefficients derivCoef <- rep( NA, nInt ) derivCoef[1] <- -2 * weights[1] * xBound[2] / ( xBound[3] - xBound[1] ) derivCoef[nInt] <- 2 * weights[nInt-1] * xBound[nInt] / ( xBound[nInt+1] - xBound[nInt-1] ) if( nInt > 2 ) { for( n in 2:( nInt-1 ) ) { derivCoef[n] <- 2 * weights[n-1] * xBound[n] / ( xBound[n+1] - xBound[n-1] ) - 2 * weights[n] * xBound[n+1] / ( xBound[n+2] - xBound[n] ) } } # variance-covariance matrix of the coefficiencts vcovCoef <- diag( xCoefSE^2 ) # standard error of the (average) semi-elasticity semElaSE <- drop( sqrt( t( derivCoef ) %*% vcovCoef %*% derivCoef ) ) # prepare object that will be returned result <- c( semEla = semElaAvg, stdEr = semElaSE ) return( result ) } @ with argument \code{xCoef} $= \delta = ( \delta_1 , \ldots , \delta_M )^{\top}$, argument \code{xShares} $= s = ( s_1 , \ldots , s_M )^{\top}$, argument \code{xBound} $= b = ( b_0 , \ldots , b_M )^{\top}$, and an optional argument \code{xCoefSE} $= \delta = ( \se( \delta_1 ), \ldots , \se( \delta_M ) )^{\top}$. <<>>= # Example ela3a <- linElaInt( xCoef = c( 0, 0.22, 0.05, 0.6 ), xShares = c( 0.35, 0.4, 0.12, 0.13 ), xBound = c( 0, 500, 1000, 1500, Inf ) ) ela3a # Example ela3b <- linElaInt( xCoef = c( 0, 0.22, 0.05, 0.6 ), xShares = c( 0.35, 0.4, 0.12, 0.13 ), xBound = c( 0, 500, 1000, 1500, Inf ), xCoefSE = c( 0, 0.002, 0.005, 0.001 ) ) ela3b @ %=================================================================================== \subsection{Effects of continuous variables when they change between discrete intervals} \label{sec:linEffInt} As in section~\ref{sec:LPMcont}, we assume a regression equation: \begin{equation} y = \beta_0 + \sum_{k=1}^K \beta_k x_k + u, \end{equation} in which all explanatory variables $x_1, \ldots , x_K$ are continuous and also all other variables and coefficients are defined as above. In this section, we derive the (approximate) effects of a variable $x_k$ on $y$ if this variable changes between $M \geq 2$~discrete intervals (e.g., age categories), e.g., from a `reference' interval~$l$ to an interval of interest~$m$: \begin{align} E_{k,ml} = & E[ y | b_{m-1} < x_k \leq b_m, x_{-k} ] - E[ y | b_{l-1} < x_k \leq b_l, x_{-k} ] \label{eq:linEffectStart} \\ = & \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus k} \beta_j x_j + \beta_k E[ x_k | b_{m-1} < x_k \leq b_m ] \\ & - \beta_0 - \sum_{j \in \{ 1, \ldots, K \} \setminus k} \beta_j x_j - \beta_k E[ x_k | b_{l-1} < x_k \leq b_l ] \nonumber \\[1ex] = & \beta_k \left( \bar{x}_{km} - \bar{x}_{kl} \right), \label{eq:linEffect} \end{align} where $x_{-k} = ( x_1 , \ldots , x_{k-1} , x_{k+1} , x_K )^{\top}$ is a vector of all explanatory variables except for~$x_k$, $b_0 , \ldots , b_M$ are the boundaries of the intervals of variable~$x_k$, and \begin{equation} \bar{x}_{km} \equiv E[ x_k | b_{m-1} < x_k \leq b_m ] \; \forall \; m = 1, \ldots , M \end{equation} are the expected values of variable~$x_k$ within specific intervals. If the expected values of variable~$x_k$ for specific intervals are unknown, it may be appropriate to approximate them by the mid-points of the respective interval boundaries (e.g., if the variable~$x_k$ has approximately a uniform distribution between the respective interval boundaries): \begin{equation} \bar{x}_{km} \approx \frac{ b_{m-1} + b_m }{2} \; \forall \; m = 1, \ldots , M. \label{eq:linEffectXBar} \end{equation} If the model specification additionally includes a quadratic term of the explanatory variable~$k$, e.g., $x_{k+1} = x_k^2$, equations~(\ref{eq:linEffectStart}) to~(\ref{eq:linEffect}) change to: \begin{align} E_{k,ml} = & E[ y | b_{m-1} < x_k \leq b_m, x_{-k} ] - E[ y | b_{l-1} < x_k \leq b_l, x_{-k} ]\\ = & \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus \{ k, k+1 \} } \beta_j x_j \\ & + \beta_k E[ x_k | b_{m-1} < x_k \leq b_m ] + \beta_{k+1} E[ x_k^2 | b_{m-1} < x_k \leq b_m ] \nonumber \\ & - \beta_0 - \sum_{j \in \{ 1, \ldots, K \} \setminus \{ k, k+1 \} } \beta_j x_j \nonumber \\ & - \beta_k E[ x_k | b_{l-1} < x_k \leq b_l ] - \beta_{k+1} E[ x_k^2 | b_{l-1} < x_k \leq b_l ] \nonumber \\[1ex] = & \beta_k \left( \bar{x}_{km} - \bar{x}_{kl} \right) + \beta_{k+1} \left( \overline{x^2_{km}} - \overline{x^2_{kl}} \right) \label{eq:quadEffect} \end{align} with \begin{equation} \overline{x^2_{km}} \equiv E[ x_k^2 | b_{m-1} < x_k \leq b_m ] \; \forall \; m = 1, \ldots , M. \end{equation} If $E[ x_k^2 | b_{m-1} < x_k \leq b_m ]$ is unknown, it may be appropriate to approximate it by assuming that variable~$x_k$ has approximately a uniform distribution between each pair of subsequent interval boundaries so that its probability density function between boundaries $b_{m-1}$ and $b_m$ is $1 / ( b_m - b_{m-1} )$: \begin{align} \overline{x^2_{km}} & \approx \int_{b_{m-1}}^{b_m} x_k^2 \; \frac{1}{b_m - b_{m-1}} \; d \, x_k\\ & = \left. \frac{1}{3} \; x_k^3 \; \frac{1}{b_m - b_{m-1}} \right|_{b_{m-1}}^{b_m} \\ & = \frac{1}{3} \; b_m^3 \; \frac{1}{b_m - b_{m-1}} - \frac{1}{3} \; b_{m-1}^3 \; \frac{1}{b_m - b_{m-1}} \\ & = \frac{ b_m^3 - b_{m-1}^3 }{ 3 \left( b_m - b_{m-1} \right)}. \label{eq:linEffectXSquaredBar} \end{align} An approximate standard error of the effect $E_{k,ml}$ calculated by equation~(\ref{eq:linEffect}) or~(\ref{eq:quadEffect}) can be obtained---as above---with the Delta method: \begin{align} \se( E_{k,ml} ) & = \sqrt{ \left( \frac{ \partial E_{k,ml} }{ \partial \left( \begin{array}{c} \beta_k \\ \beta_{k+1} \end{array} \right) } \right)^{\top} \var\left( \begin{array}{c} \beta_k \\ \beta_{k+1} \end{array} \right) \; \frac{ \partial E_{k,ml} }{ \partial \left( \begin{array}{c} \beta_k \\ \beta_{k+1} \end{array} \right) } }, \label{eq:linEffectSE} \end{align} where $\se( E_{k,ml} )$ indicates the (approximate) standard error of $E_{k,ml}$, $\var\left( ( \beta_k , \beta_{k+1} )^{\top} \right)$ indicates the variance-covariance matrix of the estimates of~$\beta_k$ and~$\beta_{k+1}$, the first element of the partial derivative of the effect~$E_{k,ml}$ w.r.t.\ the coefficients is: \begin{equation} \frac{\partial E_{k,ml}}{\partial \beta_k} = \bar{x}_{km} - \bar{x}_{kl} \label{eq:linEffectDeriv} \end{equation} and the second element of this partial derivative is zero if there is no quadratic term of~$x_k$ (so that the effect is calculated by equation~\ref{eq:linEffect}) and it is: \begin{equation} \frac{\partial E_{k,ml}}{\partial \beta_{k+1}} = \overline{x^2_{km}} - \overline{x^2_{kl}} \label{eq:quadEffectDeriv} \end{equation} if there is a quadratic term of~$x_k$ (so that the effect is calculated by equation~\ref{eq:quadEffect}). If the covariance between $\beta_k$ and $\beta_{k+1}$ is unknown, one ccould assume that it is zero. The following code defines a helper function that calculates $\overline{x^2_{km}}$ according to equation~(\ref{eq:linEffectXSquaredBar}): <<>>= EXSquared <- function( lowerBound, upperBound ) { result <- ( upperBound^3 - lowerBound^3 )/( 3 * ( upperBound - lowerBound ) ) return( result ) } @ with arguments \code{lowerBound} $= b_{m-1}$ and \code{upperBound} $= b_m$ the lower and upper boundaries of the interval, respectively. Using helper functions \code{EXSquared} (equation~\ref{eq:linEffectXSquaredBar}) and \code{elaIntBounds} (checking arguments \code{refBound} and \code{intBound}), the following function calculates the effect $E_{k,ml}$ and its approximate standard error $\se(E)$ according to equations~(\ref{eq:linEffect}), (\ref{eq:linEffectXBar}), (\ref{eq:quadEffect}), and (\ref{eq:linEffectSE}) to (\ref{eq:quadEffectDeriv}): << >>= linEffInt <- function( xCoef, refBound, intBound, xCoefSE = rep( NA, length( xCoef ) ) ){ if( ! length( xCoef ) %in% c( 1, 2 ) ) { stop( "argument 'xCoef' must be a scalar or vector with 2 elements" ) } refBound <- elaIntBounds( refBound, 1, argName = "refBound" ) intBound <- elaIntBounds( intBound, 1, argName = "intBound" ) if( length( xCoef ) != length( xCoefSE ) ) { stop( "arguments 'xCoef' and 'xCoefSE' must have the same length" ) } if( length( xCoef ) == 1 ) { xCoef <- c( xCoef, 0 ) xCoefSE <- c( xCoefSE, 0 ) } # difference between the xBars of the two intervals xDiff <- mean( intBound ) - mean( refBound ) # difference between the xSquareBars of the two intervals xSquaredDiff <- EXSquared( intBound[1], intBound[2] ) - EXSquared( refBound[1], refBound[2] ) # effect E_{k,ml} eff <- xCoef[1] * xDiff + xCoef[2] * xSquaredDiff # partial derivative of E_{k,ml} w.r.t. the beta_k and beta_{k+1} derivCoef <- c( xDiff, ifelse( xCoef[2] == 0, 0, xSquaredDiff ) ) # variance covariance of the coefficients (covariances set to zero) vcovCoef <- diag( xCoefSE^2 ) # approximate standard error of the effect effSE <- drop( sqrt( t( derivCoef ) %*% vcovCoef %*% derivCoef ) ) # object to be returned result <- c( effect = eff, stdEr = effSE ) return( result ) } @ with argument \code{xCoef} $= \beta_k$ or $( \beta_k, \beta_{k+1})^\top$ the coefficient(s) or marginal effect(s) of interest, \code{refBound} $= ( b_{l-1}, b_l )$ the boundaries of the reference interval, \code{\mbox{intBound}} $= ( b_{m-1}, b_m )$ the boundaries of the interval of interest, and optional argument \code{xCoefSE} $= \se( \beta_k )$ or $(\se( \beta_k ), \se(\beta_{k+1}))^\top$ the standard error(s) of the coefficient(s) or marginal effect(s) of interest. << >>= # Example # model equation that is linear in x_k eff1a <- linEffInt( 0.4, refBound = c( 19, 34 ), intBound = c( 35, 55 ) ) eff1a eff1b <- linEffInt( 0.4, refBound = c( 19, 34 ), intBound = c( 35, 55 ), xCoefSE = 0.03 ) eff1b # Example # model equation that is quadratic in x_k eff2a <- linEffInt( c( 0.4, -0.0003 ), refBound = c( 19, 34 ), intBound = c( 35, 55 ) ) eff2a eff2b <- linEffInt( c( 0.4, -0.0003 ), refBound = c( 19, 34 ), intBound = c( 35, 55 ), xCoefSE = c( 0.002, 0.000001 ) ) eff2b @ %===================================================================================== \subsection{Grouping and re-basing categorical variables} \label{sec:lingroup} We consider a regression model: \begin{align} y & = \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m \in \{ 1, \ldots , M \} \setminus m^* } \delta_m D_m + u,\\ D_m & = \begin{cases} 1 & \text{ if } x_k \in c_m \\ 0 & \text{ otherwise} \end{cases} \quad \forall \; m = 1, \ldots , M, \end{align} where the variable of interest, $x_k$, is a categorical variable with $M$ mutually exclusive categories $c_1, \ldots, c_M$ with $c_m \cap c_l = \emptyset \; \forall \; m \neq l$, category $c_{m^*}$ is used as reference category, and all other variables and coefficients are defined as above. For notational simplification of the following derivations, we define the (non-estimated) coefficient of the reference category to be zero, i.e., $\delta_{m^*} \equiv 0$. We want to obtain the effect of a change of variable~$x_k$ from a reference category~$c_r^*$ to a category of interest~$c_l^*$: \begin{align} E_{k,lr} & = E[ y | x_k \in c_l^*] - E[ y | x_k \in c_r^* ], \label{eq:effCat} \end{align} where categories~$c_r^*$ and/or~$c_l^*$ may comprise multiple original categories $c_1, \ldots, c_M$. Vectors $v_r = ( v_{r1} , \ldots , v_{rM} )^{\top}$ and $v_l = ( v_{l1} , \ldots , v_{lM} )^{\top}$ indicate, which of the original categories $c_1, \ldots, c_M$ are included in categories~$c_r^*$ and~$c_l^*$, respectively: \begin{align} v_{nm} & = \begin{cases} 1 & \text{ if } c_m \in c_n^* \\ 0 & \text{ if } c_m \notin c_n^* \end{cases} \forall \; m = 1, \ldots, M ; n \in \{ r , l \} \end{align} In the following, we derive the effect of a change of variable~$x_k$ from a reference category~$c_r^*$ to a category of interest~$c_l^*$, $E_{k,lr}$ as defined in equation~(\ref{eq:effCat}): \begin{align} E_{k,lr} = & E[ y | x_k \in c_l^* ] - E[ y | x_k \in c_r^* ] \\ = & \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m E[ D_m | x_k \in c_l^* ] \\ & - \beta_0 - \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j - \sum_{m=1}^{M} \delta_m E[ D_m | x_k \in c_r^* ] \nonumber \\ = & \sum_{m=1}^{M} \delta_m \left( E[ D_m | x_k \in c_l^* ] - E[ D_m | x_k \in c_r^* ] \right) \\ = & \sum_{m=1}^{M} \delta_m \left( D_{ml} - D_{mr} \right) \label{eq:linEffGr} \end{align} with \begin{align} D_{mn} & \equiv E[ D_m | x_k \in c_n^* ] \\ & = \frac{ P( D_m = 1 \cap x_k \in c_n^* ) }{ P( x_k \in c_n^* ) } \\ & = \frac{ E[ D_m ] \; v_{nm}}{ P( x_k \in c_n^* ) } \\ & = \frac{ s_m v_{nm} }{ \sum_{k=1}^M s_k v_{nk} } \label{eq:linEffGrD}\\ & \qquad \forall \; m = 1, \ldots, M ; n \in \{ r , l \}, \nonumber \end{align} where $s_m = E[ D_m ] \; \forall \; m = 1, \ldots , M$ are the shares of observations, where variable~$x_k$ is in category~$c_m$. The delta method can be used to calculate approximate standard errors of~$E_{k,lr}$. When setting all covariances between the coefficients to zero, an approximate standard error of~$E_{k,lr}$ can be obtained by: \begin{align} \se\left( E_{k,lr} \right) \approx \sqrt{ \sum_{m=1}^{M} \left( \se\left( \delta_m \right) \; \left( D_{ml} - D_{mr} \right) \right)^2 }, \label{eq:linEffGrSE} \end{align} where $\se\left( \delta_m \right)$ are the standard errors of the estimates of $\delta_m$. The standard error of the non-estimated coefficient of the reference category~$m^*$ is set to zero, i.e., $\se\left( \delta_{m^*} \right) = 0$. As an illustrative example, we assume that variable~$x_k$ has $M=5$ categories, whereas the first category is used as reference category in the regression model, i.e., $m^* = 1$ and $\delta_1 \equiv 0$, and we want to measure the effect of variable~$x_k$ changing from the combined reference category~$\{3,4\}$ to the combined category of interest~$\{1,2\}$ so that $v_r = ( 0, 0, 1, 1, 0 )^{\top}$ and $v_l = ( 1, 1, 0, 0, 0 )^{\top}$: \begin{align} E_{k,\{1,2\}\,\{3,4\}} = & \delta_1 \left( \frac{ s_1 }{ s_1 + s_2 } - \frac{ 0 }{ s_3 + s_4 } \right) + \delta_2 \left( \frac{ s_2 }{ s_1 + s_2 } - \frac{ 0 }{ s_3 + s_4 } \right)\\ & + \delta_3 \left( \frac{ 0 }{ s_1 + s_2 } - \frac{ s_3 }{ s_3 + s_4 } \right) + \delta_4 \left( \frac{ 0 }{ s_1 + s_2 } - \frac{ s_4 }{ s_3 + s_4 } \right) \nonumber\\ & + \delta_5 \left( \frac{ 0 }{ s_1 + s_2 } - \frac{ 0 }{ s_3 + s_4 } \right) \nonumber\\ = & \frac{ s_2 \cdot \delta_2 }{ s_1 + s_2 } - \frac{ s_3 \cdot \delta_3 }{ s_3 + s_4 } - \frac{ s_4 \cdot \delta_4 }{ s_3 + s_4 }. \label{eq:effectRB} \end{align} with approximate standard error: \begin{align} \se\left( E_{k,\{1,2\}\,\{3,4\}} \right) &= \sqrt{ \left( \frac{ s_2 \cdot \se(\delta_2) }{ s_1 + s_2 } \right)^2 - \left( \frac{ s_3 \cdot \se(\delta_3) }{ s_3 + s_4 } \right)^2 - \left( \frac{ s_4 \cdot \se(\delta_4) }{ s_3 + s_4 } \right)^2 } . \label{eq:effectRBSE} \end{align} The following function calculates the effect~$E_{k,lr}$ and the its standard error according to equations~(\ref{eq:linEffGr}), (\ref{eq:linEffGrD}), and~(\ref{eq:linEffGrSE}): << >>= linEffGr <- function( xCoef, xShares, Group, xCoefSE = rep( NA, length( xCoef ) ) ){ if( sum( xShares ) > 1 ){ stop( "the shares in argument 'xShares' sum up to a value larger than 1" ) } if( length( xCoef ) != length( xShares ) ){ stop( "arguments 'xCoef' and 'xShares' must have the same length" ) } if( length( xCoef ) != length( Group ) ){ stop( "arguments 'xCoef' and 'Group' must have the same length" ) } if( ! all( Group %in% c( -1, 0, 1 ) ) ){ stop( "all elements of argument 'Group' must be -1, 0, or 1" ) } # D_mr DRef <- xShares * ( Group == -1 ) / sum( xShares[ Group == -1 ] ) # D_ml DEffect <- xShares * ( Group == 1 ) / sum( xShares[ Group == 1 ] ) # effect: sum of delta_m * ( D_ml - D_mr ) effeG <- sum( xCoef * ( DEffect - DRef ) ) # approximate standard error effeGSE <- sqrt( sum( ( xCoefSE * ( DEffect - DRef ) )^2 ) ) result <- c( effect = effeG, stdEr = effeGSE ) return( result ) } @ with argument \code{xCoef} $= (\delta_1, \ldots, \delta_{M})^{\top}$ a vector of coefficients of the categorical variable of interest, \textbf{where the coefficient of the reference group is set to 0}, argument \code{xShares} $= ( s_1, \ldots, s_M )^{\top}$ a vector of the corresponding shares of each category, argument \code{Group} $= ( D_{1l} - D_{1r}, \ldots , D_{Ml} - D_{Mr} )^{\top}$ a vector consisting of $-1$s, 0s, and 1s, where a $-1$ indicates that the category belongs to the reference category, a 1 indicates that the category belongs to the category of interest, and a zero indicates that the category is neither in the base category nor in the category of interest, and optional argument \code{xCoefSE} $= ( \se(\delta_1), \ldots, \se(\delta_M))^{\top}$, \textbf{where the standard error of the coefficient of the reference group is set to 0}. Elements of arguments \code{xCoef}, \code{xShares}, \code{Group}, and \code{xCoefSE} that belong to a category that is neither in the reference category nor in the category of interest, i.e., categories~$m$ with $D_{mr} = D_{ml} = 0$, can be omitted but the ommission must be consistent for all four arguments. << >>= # Example: eff3a <- linEffGr( xCoef = c( 0, 0.2, 0.04, 0.06, 0.3 ), xShares = c( 0.14, 0.35, 0.3, 0.01, 0.2 ), Group = c( -1, 1, -1, -1, 0 ) ) eff3a # Example: eff3b <- linEffGr( xCoef = c( 0, 0.2, 0.04, 0.06, 0.3 ), xShares = c( 0.14, 0.35, 0.3, 0.01, 0.2 ), Group = c( -1, 1, -1, -1, 0 ), xCoefSE = c( 0, 0.0001, 0.002, 0.05, 0.09 )) eff3b # Example: eff3c <- linEffGr( xCoef = c( 0, 0.2, 0.04, 0.06 ), xShares = c( 0.14, 0.35, 0.3, 0.01 ), Group = c( -1, 1, -1, -1 )) eff3c # Example: eff3d <- linEffGr( xCoef = c( 0, 0.2, 0.04, 0.06 ), xShares = c( 0.14, 0.35, 0.3, 0.01 ), Group = c( -1, 1, -1, -1 ), xCoefSE = c( 0, 0.0001, 0.002, 0.05 )) eff3d @ %============================================================================================== \section{Binary probit model, multivariate probit model, and ordered probit model} \subsection{Semi-elasticities from continuous explanatory variables (linear and quadratic)} \label{sec:probEla} The (binary) probit model with continuous explanatory variables is specified as: \begin{equation} \Pr( y = 1 | x ) = \Phi\left( \beta_0 + \sum_{j=1}^K \beta_j x_j \right) , \label{eq:probit} \end{equation} where $y \in \{0, 1 \}$ is a binary dependent variable, $x = (x_1, \ldots, x_K )^{\top}$ is a vector of $K$ continuous explanatory variables, $\beta = ( \beta_0, \ldots, \beta_K )^{\top}$ is a vector of $K + 1$ unknown coefficients, and $\Phi(\cdot)$ is the cumulative distribution function of the standard normal distribution. The semi-elasticity of the $k$th (continuous) explanatory variable is: \begin{equation} \epsilon_k = \frac{\partial \Pr( y = 1 )}{\partial x_k } \cdot x_k = \phi \left( \beta_0 + \sum_{j = 1}^K \beta_j x_j \right) \cdot \beta_k x_k, \label{eq:probEla} \end{equation} where $\phi(\cdot)$ is the probability density function of the standard normal distribution. The interpretation of the semi-elasticity~$\epsilon_k$ is identical to the one described in section~\ref{sec:LPMcont}. The value of the semi-elasticity~$\epsilon_k$ depends on the values of all explanatory variables~$x_1 , \ldots , x_K$. In order to obtain a `representative' single value of this semi-elasticity, it is usually calculated at the sample mean values of the explanatory variables $\bar{x}_1 , \ldots , \bar{x}_K$. In case of binomial or multinomial probit models with two or more dependent variables $y_1 , \ldots , y_N$, we only calculate the marginal effects of $x_k$ on the unconditional expectations $E[y_n | x_1 , \ldots , x_K ]$ (and disregard marginal effects on the conditional expectations $E[ y_n | x_1 , \ldots , x_K, y_1 , \ldots , y_{n-1} , y_{n+1} , \ldots , y_K$). This has the advantage that we can ignore the interrelations between the regression models for the different dependent variables and obtain the semi-elasticities in the same way as for the univariate probit model, i.e., we can ignore the coefficients of the regression models for the other dependent variables as well as the coefficients of correlation between the error terms. Hence, equation~(\ref{eq:probEla}) still applies. The ordered probit model, where the dependent variable can have $P$~distinct and strictly ordered values, i.e., $y \in \{ 1, \ldots , P \}$, can be specified as: \begin{equation} \Pr( y = p ) = \Phi\left( \mu_p - \sum_{j = 1}^K \beta_j x_j \right) - \Phi\left( \mu_{p - 1} - \sum_{j = 1}^K \beta_j x_j \right) \; \forall \; p \in \{ 1, \ldots , P \}, \end{equation} where $\mu_0 , \ldots , \mu_P$ are the break points, of which $\mu_0 = - \infty$, $\mu_P = \infty$, and $\mu_1 , \ldots , \mu_{P-1}$ are unknown and, thus, need to be estimated. We create a new binary dependent variable~$y^*$ by dividing the $P$~distinct values of the dependent variable~$y$ into two categories: \begin{equation} y^* = \begin{cases} 0 & \text{ if } y \in \{ 1, \ldots , p^* - 1 \} \\ 1 & \text{ if } y \in \{ p^* , \ldots , P \} \end{cases} \end{equation} with $p^* \in \{ 2, \ldots, P \}$. This reduces the ordered probit model to a binary probit model: \begin{align} \Pr( y^* = 1 ) & = \Pr( y \in \{ p^*, \ldots , P \} ) \\ & = \sum_{p = p^*}^{P} \Pr( y = p ) \\ & = \sum_{p = p^*}^{P} \left( \Phi\left( \mu_p - \sum_{j = 1}^K \beta_j x_j \right) - \Phi\left( \mu_{p - 1} - \sum_{j = 1}^K \beta_j x_j \right) \right)\\ & = \sum_{p = p^*}^{P} \Phi\left( \mu_p - \sum_{j = 1}^K \beta_j x_j \right) - \sum_{p = p^*}^{P} \Phi\left( \mu_{P - 1} - \sum_{j = 1}^K \beta_j x_j \right) \\ & = \sum_{p = p^*}^{P} \Phi\left( \mu_p - \sum_{j = 1}^K \beta_j x_j \right) - \sum_{p = p^* - 1}^{P-1} \Phi\left( \mu_{p} - \sum_{j = 1}^K \beta_j x_j \right) \\ & = \Phi\left( \mu_P - \sum_{j = 1}^K \beta_j x_j \right) - \Phi\left( \mu_{p^* - 1} - \sum_{j = 1}^K \beta_j x_j \right) \\ & = \Phi\left( \infty \right) - \Phi\left( \mu_{p^* - 1} - \sum_{j = 1}^K \beta_j x_j \right) \\ & = 1 - \Phi\left( \mu_{p^* - 1} - \sum_{j = 1}^K \beta_j x_j \right) \\ & = \Phi\left( - \mu_{p^* - 1} + \sum_{j = 1}^K \beta_j x_j \right), \end{align} where the intercept of the binary probit model is equal to the negative value of the break point that separates $y^* = 0$ from $y^* = 1$, i.e., $\beta_0 = -\mu_{p^* - 1}$.% \footnote{ If the ordered probit model is estimated with intercept, say, $\beta_0^*$ and (for identification) by normalising the first (internal) break point to zero, i.e., $\mu_1 = 0$, the ordered probit model can be simplified to a binary probit model with intercept $\beta_0 = \beta_0^* - \mu_{p^* - 1}$. } Hence, we can derive the semi-elasticity in the same way as for the binary probit model, i.e., by applying equation~(\ref{eq:probEla}). If the model specification additionally includes a quadratic term of the explanatory variable $k$, e.g., $x_{k+1} = x^2_k$, the semi-elasticity of this explanatory variable is: \begin{equation} \epsilon_k = \frac{\partial \Pr( y = 1 )}{\partial x_k} \cdot x_k = \phi \left( \beta_0 + \sum_{j = 1}^K \beta_j x_j \right) \cdot ( \beta_k x_k + 2 \beta_{k+1} x_k^2 ), \label{eq:probEla2} \end{equation} In order to calculate the standard errors of the marginal effects (and semi-elasticities) of a probit model, it is common to apply the Delta method as indicated in equation~(\ref{eq:linElaQuadSE}). But as in most cases the variance-covariance matrix of the regression coefficients is unknown, standard errors can only be approximated through a simplified Delta method that sets all covariances to zero. However, we conducted simulations that show that setting all covariances to zero leads to standard errors that are much larger than correctly calculated standard errors (i.e.\ using the estimated covariances rather than setting them to zero). In contrast, our simulations indicated that simplifying the calculation of the standard errors by assuming that $\phi \left( \beta_0 + \sum_{j = 1}^K \beta_j x_j \right)$ is a constant, i.e., assuming $\partial \phi \left( \beta_0 + \sum_{j = 1}^K \beta_j x_j \right) / \partial \beta_n = 0 \; \forall \; n = 0, \ldots , K$, gives standard errors that are much closer to correctly calculated standard errors than including the full gradient vector of equation~(\ref{eq:probEla}) with respect to the coefficients. Under this simplifying assumption, only one element of the gradient of equation~(\ref{eq:probEla}) with respect to the coefficients is non-zero: \begin{align} \frac{\partial \epsilon_k}{\partial \beta_k} &= \phi \left( \beta_0 + \sum_{j = 1}^K \beta_j x_j \right) \cdot x_k \end{align} and an approximate standard error of $\epsilon_k$ can be obtained by: \begin{align} \se( \epsilon_k ) &= \sqrt{ \phi \left( \beta_0 + \sum_{j = 1}^K \beta_j x_j \right) x_k \cdot \var( \beta_k ) \cdot \phi \left( \beta_0 + \sum_{j = 1}^K \beta_j x_j \right) x_k } \label{eq:probse1}\\ &= \phi \left( \beta_0 + \sum_{j = 1}^K \beta_j x_j \right) \cdot x_k \cdot \se( \beta_k ), \label{eq:probse} \end{align} where $\se( \beta_k )$ is the standard error of the coefficient of interest. If the variable of interest enters the regression equation in quadratic form as defined in equation~(\ref{eq:probEla2}), the simplified gradient vector has two non-zero elements: \begin{align} \frac{\partial \epsilon_k} { \partial \begin{pmatrix} \beta_k \\ \beta_{k + 1 } \end{pmatrix}} &= \begin{pmatrix} \phi \left( \beta_0 + \sum_{j = 1}^K \beta_j x_j \right) \cdot x_k \\ \\ \phi \left( \beta_0 + \sum_{j = 1}^K \beta_j x_j \right) \cdot 2 x_k^2 \end{pmatrix} \label{eq:probGrad2} \end{align} and an approximate standard error can be calculated by: \begin{align} \se( \epsilon_k ) &= \sqrt{ \left( \frac{ \partial \epsilon_k } { \partial \begin{pmatrix} \beta_k \\ \beta_{k + 1 } \end{pmatrix}} \right)^{\top} \cdot \var \begin{pmatrix} \beta_k & 0 \\ 0 & \beta_{k+1} \end{pmatrix} \cdot \left( \frac{ \partial \epsilon_k } { \partial \begin{pmatrix} \beta_k \\ \beta_{k + 1 } \end{pmatrix}} \right) } . \label{eq:probse2} \end{align} % helper function (code not shown here but in the appendix) <>= checkXPos <- function( xPos, minLength, maxLength, minVal, maxVal, requiredVal = NA ) { if( any( xPos != round( xPos ) ) ) { stop( "argument 'xPos' must be a vector of integers" ) } if( length( xPos ) < minLength ) { stop( "argument 'xPos' must have a length equal to or larger than ", minLength ) } if( length( xPos ) > maxLength ) { stop( "argument 'xPos' must have a length smaller than or equal to ", maxLength ) } if( any( xPos < minVal ) ) { stop( "all elements of argument 'xPos' must be equal to or larger than ", minVal ) } if( any( xPos > maxVal ) ) { stop( "all elements of argument 'xPos' must be smaller than or equal to ", maxVal ) } if( max( table( xPos ) ) > 1 ) { stop( "all elements of argument 'xPos' may only occur once" ) } if( !is.na( requiredVal ) ) { if( sum( xPos == requiredVal ) != 1 ) { stop( "argument 'xPos' must have exactly one element that is ", requiredVal ) } } } @ % helper function (code not shown here but in the appendix) <>= checkXBeta <- function( xBeta ) { if( any( abs( xBeta ) > 3.5 ) ) { warning( "At least one x'beta has an implausible value: ", paste( xBeta, collapse = ", " ) ) } } @ Using the helper functions \code{checkXPos} (checking argument \code{xPos}, see appendix~\ref{sec:helperFunctions}) and \code{checkXBeta} (checking $\beta_0 + \sum_{j=1}^K \beta_j x_j$ for plausible values, see appendix~\ref{sec:helperFunctions}), the following function calculates the semi-elasticity, $\epsilon_k$, and the corresponding standard error according to equations~(\ref{eq:probEla}), (\ref{eq:probEla2}), (\ref{eq:probGrad2}), and (\ref{eq:probse2}): << >>= probEla <- function( allCoef, allXVal, allCoefSE = rep( NA, length( allCoef )), xPos ){ nCoef <- length( allCoef ) if( nCoef != length( allCoefSE ) ) { stop( "arguments 'allCoef' and 'allCoefSE' must have the same length" ) } if( length( allCoef ) != length( allXVal ) ) { stop( "arguments 'allCoef' and 'allXVal' must have the same length" ) } checkXPos( xPos, minLength = 1, maxLength = 2, minVal = 1, maxVal = nCoef ) if( length( xPos ) == 2 ){ xCoef <- allCoef[ xPos ] xCoefSE <- allCoefSE[ xPos ] if( !isTRUE( all.equal( allXVal[ xPos[2] ], allXVal[ xPos[1] ]^2 ) ) ) { stop( "the value of 'allXVal[ xPos[2] ]' must be equal", "to the squared value of 'allXVal[ xPos[1] ]' " ) } } else if( length( xPos ) == 1 ) { xCoef <- c( allCoef[ xPos ], 0 ) xCoefSE <- c( allCoefSE[ xPos ], 0 ) } else { stop( "argument 'xPos' must be a scalar or a vector with two elements" ) } xVal <- allXVal[ xPos[ 1 ] ] xBeta <- sum( allCoef * allXVal ) checkXBeta( xBeta ) dfun <- pnorm( xBeta ) semEla <- ( xCoef[ 1 ] + 2 * xCoef[ 2 ] * xVal ) * xVal * dfun derivCoef <- c( dfun * xVal, ifelse( length( xPos ) == 1, 0, dfun * 2 * xVal^2 ) ) vcovCoef <- diag( xCoefSE^2 ) semElaSE <- drop( sqrt( t( derivCoef ) %*% vcovCoef %*% derivCoef ) ) result <- c( semEla = semEla, stdEr = semElaSE ) return( result ) } @ with argument \code{allCoef} $= ( \beta_0 , \ldots , \beta_K )^{\top}$, or alternatively in the case of an ordered probit model \code{allCoef} $= (-\mu_{p^* - 1}, \beta_1, \ldots, \beta_K)^\top$, a vector of all coefficients; argument \code{allXVal} $= ( 1, x_1 , \ldots , x_K )^{\top}$ a vector of the values of all corresponding explanatory variables including the intercept (e.g., their sample means); argument \code{xPos} = $k+1$ or $( k+1 , k+2 )^{\top}$ a scalar or vector with two elements that indicates the position of the variable of interest and its coefficient and eventually the position of the squared variable of interest and its coefficient in vectors \code{allXVal} and \code{allCoef};% \footnote{ If argument \code{xPos} is a scalar, equation~(\ref{eq:probEla2}) simplifies to equation~(\ref{eq:probEla}) and equations~(\ref{eq:probGrad2}) and~(\ref{eq:probse2}) simplify to equation~(\ref{eq:probse}). } and optional argument \code{allCoefSE} $= ( \se( \beta_0 ), \ldots , \se( \beta_K ) )^\top$ or alternatively in the case of an ordered probit model \code{allCoefSE} $= ( \se( \mu_{m^* - 1}), \se( \beta_1 ), \ldots, \se( \beta_K ))^\top$,% \footnote{ Note that $\se( - \mu_{m^* - 1}) = \sqrt{ \var( - \mu_{m^* - 1} ) } = \sqrt{ ( \partial ( - \mu_{m^* - 1} ) / \partial \mu_{m^* - 1} )^2 \; \var( \mu_{m^* - 1} ) } = \sqrt{ ( -1 )^2 \var( \mu_{m^* - 1} ) } = \sqrt{ \var( \mu_{m^* - 1} ) } = \se( \mu_{m^* - 1})$. } the standard errors of the coefficients. << >>= # Example ela4a <- probEla( c( 0.445, 0.03, 0.00002, 0.067, 0.89, 0.124 ), c( 1, 2.34, 3.3, 3.3^2, 0.0, 0.987 ), xPos = 2 ) ela4a # Example ela4b <- probEla( c( 0.445, 0.03, 0.00002, 0.067, 0.89, 0.124 ), c( 1, 2.34, 3.3, 3.3^2, 0.0, 0.987 ), c( 0.032, 0.004, 0.00001, 0.034, 0.0009, 0.056 ), xPos = 2 ) ela4b # Example ela4c <- probEla( c( 0.445, 0.03, 0.00002, 0.067, 0.89, 0.124 ), c( 1, 2.34, 3.3, 3.3^2, 0.0, 0.987 ), c( 0.032, 0.004, 0.00001, 0.034, 0.0009, 0.056 ), xPos = c( 3, 4 )) ela4c @ %========================================================================================== \subsection{Semi-elasticities from interval-coded explanatory variables} \label{sec:probInt} Probit model, where the $k$th explanatory variable is interval-coded: \begin{align} \Pr( y = 1 | x ) &= \Phi\left( \beta_0 + \sum_{j\in\{1, \ldots, K \} \setminus k} \beta_j x_j + \sum_{m \in \{1,...,M \} \setminus m^*} \delta_m D_m \right) \label{eq:probInterv}\\ D_m &= \begin{cases} 1 & \text{ if } b_{m-1} < x_k \leq b_m \\ 0 & \text{ otherwise} \end{cases} \quad \forall \; m = 1, \ldots , M, \end{align} where all variables and coefficients are as described in section~\ref{sec:linInterval}. As described in section~\ref{sec:linInterval}, the semi-elasticities of the interval-coded variable~$x_k$ can be calculated as weighted averages by equation~(\ref{eq:linElaIntAvg}): \begin{align} \epsilon_k \equiv \frac{ \partial E[y] }{ \partial x_k } \; x_k \approx \sum_{m=1}^{M-1} w_m \; \epsilon_{km} \label{eq:probElaIntAvg} \end{align} with weights~$w_1, \ldots, w_{M-1}$ as defined in equation~(\ref{eq:linElaIntWeights}). However, in contrast to section~\ref{sec:linInterval}, the numerators of the semi-elasticities~$\epsilon_{km}$ around each inner boundary~$b_1, \ldots , b_{M-1}$ defined in equation~(\ref{eq:linElaIntBound}) cannot be simplified to differences between coefficients so that the derivations from equation~(\ref{eq:linElaIntBound}) to equation~(\ref{eq:linElaInt}) become: \begin{align} \epsilon_{km} \approx \; & \frac{ E[ y | b_m < x_k \leq b_{m+1} ] - E[ y | b_{m-1} < x_k \leq b_m ] }{ E[ x_k | b_m < x_k \leq b_{m+1} ] - E[ x_k | b_{m-1} < x_k \leq b_m ] } \; b_m \\[1ex] \approx \; & 2 \; \frac{ \Phi \left( \beta_0 + \sum_{j\in\{1, \ldots, K \} \setminus k} \beta_j x_j + \delta_{m+1} \right) -\Phi \left( \beta_0 + \sum_{j\in\{1, \ldots, K \} \setminus k} \beta_j x_j + \delta_m \right) }{ b_{ m+1 } - b_{ m-1 }} \cdot b_m \label{eq:probElaInt}\\ & \forall \; m = 1, \ldots, M-1 . \nonumber \end{align} An approximate standard error of the semi-elasticity of interval-coded explanatory variables of probit models defined in equations~(\ref{eq:probElaIntAvg}) and~(\ref{eq:probElaInt}) can be obtained by using the Delta-method: \begin{align} \se( \epsilon_k ) & = \sqrt{ \frac{ \partial \epsilon_k }{ \partial \left( \beta^{\top} \; \delta^{\top} \right) } \; \var \left( \begin{array}{c} \beta \\ \delta \end{array} \right) \; \frac{ \partial \epsilon_k }{ \partial \left( \begin{array}{c} \beta \\ \delta \end{array} \right) } }, \label{eq:probElaIntSE} \end{align} where $\se( \epsilon_k )$ indicates the (approximate) standard error of $\epsilon_k$, $\var \left( \left( \beta^{\top} \; \delta^{\top} \right)^{\top} \right)$ indicates the variance-covariance matrix of the estimates of the coefficient vector~% $\left( \beta^{\top} \; \delta^{\top} \right)^{\top}$, and the gradient vector is: \begin{equation} \frac{ \partial \epsilon_k} { \partial \begin{pmatrix} \beta \\ \delta \end{pmatrix}} = \sum_{m = 1}^{M-1} w_m \frac{ \partial \epsilon_{km}} { \partial \begin{pmatrix} \beta \\ \delta \end{pmatrix}} \end{equation} with the elements of the gradient vectors~% $\partial \epsilon_{km} / \partial ( \beta^{\top} \delta^{\top} )^{\top}$: \begin{align} \frac{ \partial \epsilon_{km} }{ \partial \beta_0 } & = \frac{2 \; ( \phi_{m+1} - \phi_m ) \; b_m}{b_{m+1} - b_{m-1}} \label{eq:probInt} \\ \frac{ \partial \epsilon_{km} }{ \partial \beta_j } & = \frac{2 \; ( \phi_{m+1} - \phi_m ) \; x_j \; b_m}{b_{m+1} - b_{m-1}} \; \forall \; j \in \{ 1, \ldots, K \} \setminus k \\ \frac{ \partial \epsilon_{km} }{ \partial \delta_m } & = \frac{- 2 \; \phi_{m} \; b_m}{b_{m+1} - b_{m-1}} \\ \frac{ \partial \epsilon_{km} }{ \partial \delta_{m+1} } & = \frac{2 \; \phi_{m+1} \; b_m}{b_{m+1} - b_{m-1}} \\ \frac{ \partial \epsilon_{km} }{ \partial \delta_n } & = 0 \; \forall \; n \in \{ 1, \ldots , M \} \setminus \{ m, m+1 \} \end{align} with \begin{align} \phi_m & = \phi \left( \beta_0 + \sum_{ j \in \{1, \ldots, K \} \setminus k } \beta_j x_j + \delta_m \right) \; \forall \; m = 1, \ldots, M \label{eq:probElaIntPhiM} \end{align} In the case of the ordered probit, one can simply replace $\beta_0$ by $-\mu_{p^* - 1}$ in equation~(\ref{eq:probElaIntPhiM}) (see section~\ref{sec:probEla}): \begin{align} \phi_m & = \phi \left( - \mu_{p^* - 1 } + \sum_{ j \in \{1, \ldots, K \} \setminus k } \beta_j x_j + \delta_m \right) \; \forall \; m = 1, \ldots, M. \end{align} Using the helper functions \code{elaIntWeights} (equation~\ref{eq:linElaIntWeights}), \code{elaIntBounds} (equation~\ref{eq:linElaIntBM}), \code{checkXPos} (checking argument \code{xPos}, see appendix~\ref{sec:helperFunctions}), and \code{checkXBeta} (checking $\beta_0 + \sum_{j\in\{1, \ldots, K \} \setminus k} \beta_j x_j + \sum_{m \in \{1,...,M \} \setminus m^*} \delta_m D_m$ for plausible values, see appendix~\ref{sec:helperFunctions}), and function \code{linElaInt} (section~\ref{sec:linInterval}), the following code defines a function that calculates the semi-elasticity defined in equations~(\ref{eq:probElaIntAvg}) and~(\ref{eq:probElaInt}) and its approximate standard error defined in equations~(\ref{eq:probElaIntSE}) to~(\ref{eq:probElaIntPhiM}): << >>= probElaInt <- function( allCoef, allXVal, xPos, xBound, allCoefSE = rep( NA, length( allCoef ) ) ){ # number of coefficients nCoef <- length( allCoef ) # number of intervals nInt <- length( xPos ) # checking arguments if( length( allXVal ) != nCoef ) { stop( "arguments 'allCoef' and 'allXVal' must have the same length" ) } if( length( allCoefSE ) != nCoef ) { stop( "arguments 'allCoef' and 'allCoefSE' must have the same length" ) } checkXPos( xPos, minLength = 2, maxLength = nCoef, minVal = 0, maxVal = nCoef, requiredVal = 0 ) if( any( allXVal[ xPos ] < 0 ) ) { stop( "all elements of argument 'allXVal'", " that are indicated by argument 'xPos'", " (i.e., the shares of observations in each interval)", " must be non-negative" ) } if( sum( allXVal[ xPos ] > 1 ) ) { stop( "the sum of the elements of argument 'allXVal'", " that are indicated by argument 'xPos'", " (i.e., the shares of observations in each interval)", " must not be larger than one" ) } # check 'xBound' and replace infinite values xBound <- elaIntBounds( xBound, nInt ) # vector of probabilities of y=1 for each interval and # vector of shares of observations in each interval xBeta <- shareVec <- rep( NA, nInt ) for( i in 1:nInt ){ allXValTemp <- replace( allXVal, xPos, 0 ) if( xPos[i] != 0 ) { allXValTemp[ xPos[i] ] <- 1 shareVec[ i ] <- allXVal[ xPos[ i ] ] } xBeta[ i ] <- sum( allCoef * allXValTemp ) } shareVec[ xPos == 0 ] <- 1 - sum( shareVec[ xPos != 0 ] ) checkXBeta( xBeta ) phiVec <- pnorm( xBeta ) # weights weights <- elaIntWeights( shareVec ) # calculation of the semi-elasticity semEla <- linElaInt( phiVec, shareVec, xBound ) ### calculation of its standard error # partial derivatives of each semi-elasticity around each boundary # w.r.t. all estimated coefficients gradM <- matrix( 0, nCoef, nInt - 1 ) gradPhiVec <- dnorm( xBeta ) for( m in 1:( nInt - 1 ) ) { gradM[ -xPos, m ] <- 2 * ( gradPhiVec[m+1] - gradPhiVec[m] ) * allXVal[ -xPos ] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) gradM[ xPos[m], m ] <- - 2 * gradPhiVec[m] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) gradM[ xPos[m+1], m ] <- 2 * gradPhiVec[m+1] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) } # partial derivative of the semi-elasticity # w.r.t. all estimated coefficients derivCoef <- rep( 0, nCoef ) for( m in 1:( nInt - 1 ) ){ derivCoef <- derivCoef + weights[m] * gradM[ , m ] } # variance-covariance matrix of the coefficiencts vcovCoef <- diag( allCoefSE^2 ) # standard error of the (average) semi-elasticity semElaSE <- drop( sqrt( t( derivCoef ) %*% vcovCoef %*% derivCoef ) ) # prepare object that will be returned result <- c( semEla[1], stdEr = semElaSE ) return( result ) } @ where argument \code{allCoef} $= ( \beta_0, \ldots, \beta_{k-1},$ $\delta_1, \ldots, \delta_{m^*-1},$ $\delta_{m^*+1}, \ldots, \delta_M,$ $\beta_{k+1}, \ldots, \beta_K )^{\top}$, or alternatively in the case of the ordered probit model \code{allCoef} $= ( -\mu_{p^* - 1}, \beta_1, \ldots, \beta_{k-1},$ $\delta_1, \ldots, \delta_{m^*-1},$ $\delta_{m^*+1}, \ldots, \delta_M,$ $\beta_{k+1}, \ldots, \beta_K )^{\top}$, specifies the vector of all coefficients; argument \code{allXVal} $= ( 1, x_1, \ldots, x_{k-1},$ $s_1, \ldots, s_{m^*-1}, s_{m^*+1}, \ldots, s_M,$ $x_{k+1}, \ldots, x_K )^{\top}$ is a vector of corresponding values of the explanatory variables (including shares of observations in each interval of variable~$x_k$ except for the reference interval~$m^*$); argument \code{xPos} $= ( k+1, \ldots, k + m^* - 1, 0, k + m^* , \ldots , k + M - 1 )$ is a vector indicating the positions of the coefficients~$\delta_1, \ldots, \delta_M$ and shares~$s_1, \ldots, s_M$ of each interval in the vectors \code{allCoef} and \code{allXVal}, whereas the position of the reference interval~$m^*$ is set to zero; argument \code{xBound} $= ( b_0 , \ldots , b_M )^{\top}$ indicates the boundaries of the intervals as in function \code{linElaInt}; and optional argument \code{allCoefSE} $= ( \se( \beta_0 ), \ldots, \se( \beta_{k-1} ),$ $\se( \delta_1 ), \ldots, \se( \delta_{m^*-1} ),$ $\se( \delta_{m^*+1} ), \ldots, \se( \delta_M ),$ $\se( \beta_{k+1} ), \ldots, \se( \beta_K ) )^{\top}$, or alternatively in the case of an ordered probit model \code{allCoefSE} $= ( \se( \mu_{m^* - 1}), \se( \beta_1 ), \ldots, \se( \beta_{k-1} ),$ $\se( \delta_1 ), \ldots, \se( \delta_{m^*-1} ),$ $\se( \delta_{m^*+1} ), \ldots, \se( \delta_M ),$ $\se( \beta_{k+1} ), \ldots, \se( \beta_K ) )^{\top}$, specifies the standard errors of all coefficients in \code{allCoef}. << >>= # Example ela5a <- probElaInt( allCoef = c( 0.33, 0.22, 0.05, 0.6 ), allXVal = c( 1, 0.4, 0.12, 0.13 ), xPos = c( 2, 0, 3, 4 ), xBound = c( 0, 500, 1000, 1500, Inf )) ela5a # Example ela5b <- probElaInt( allCoef = c( 0.33, 0.22, 0.05, 0.6 ), allXVal = c( 1, 0.4, 0.12, 0.13 ), xPos = c( 2, 0, 3, 4 ), xBound = c( 0, 500, 1000, 1500, Inf ), allCoefSE = c( 0.003, 0.045, 0.007, 0.009 )) ela5b @ %========================================================================================== \subsection{Effects of continuous variables when they change between discrete intervals} \label{sec:probEffInt} As in section~\ref{sec:probEla} we assume the following model: \begin{equation} \Pr( y = 1 | x ) = \Phi\left( \beta_0 + \sum_{j=1}^K \beta_j x_j \right) , \end{equation} In order to derive the (approximate) effect of a continuous variable $x_k$ on $Pr( y=1 )$ given that $x_k$ changes between $M \geq 2$ intervals, we modify equation~(\ref{eq:linEffect}) into: \begin{align} E_{k,ml} = \; & E[ y | b_{m-1} < x_k \leq b_m ] - E[ y | b_{l-1} < x_k \leq b_l ] \label{eq:probEffStart}\\ \approx \; & \Phi\left( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_j x_j + \beta_k \; E[ x_k | b_{m-1} < x_k \leq b_m ] \right) \label{eq:probEffMiddle}\\ & - \Phi\Bigg( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_j x_j + \beta_k \; E[ x_k | b_{l-1} < x_k \leq b_l ] \Bigg) \nonumber \\ = \; & \Phi\left( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_j x_j + \beta_k \bar{x}_{km} \right) \label{eq:probEff} \\ &- \Phi\left( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_j x_j + \beta_k \bar{x}_{kl} \right), \nonumber \end{align} where $\bar{x}_{km} \equiv E[ x_k | b_{m-1} < x_k \leq b_m ] \; \forall \; m = 1, \ldots, M$ can be approximated by the mid-points of the intervals as indicated by equation~(\ref{eq:linEffectXBar}).% \footnote{\label{foot:PhiNonLin} The derivation from equation~(\ref{eq:probEffStart}) to equation~(\ref{eq:probEffMiddle}) is approximate only. This derivation would be exact if the cumulative distribution function of the normal distribution~% $\Phi(\cdot)$ would be linear but this function is clearly non-linear. } For model specifications that additionally include a quadratic term of the explanatory variable, equation~(\ref{eq:quadEffect}) modifies to: \begin{align} E_{k,ml} = \; & E[ y | b_{m-1} < x_k \leq b_m, x_{-k} ] - E[ y | b_{l-1} < x_k \leq b_l, x_{-k} ] \label{eq:probQuadEffStart}\\ \approx \; & \Phi \Bigg( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus \{ k, k+1 \}} \beta_j x_j \label{eq:probQuadEffMiddle} \\ & \hspace{6mm} + \beta_k E[ x_k | b_{m-1} < x_k \leq b_m ] + \beta_{k+1} E[ x_k^2 | b_{m-1} < x_k \leq b_m ] \Bigg) \nonumber \\ &- \Phi \Bigg( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus \{ k, k+1 \}} \beta_j x_j \nonumber \\ & \hspace{8mm} + \beta_k E[ x_k | b_{l-1} < x_k \leq b_l ] + \beta_{k+1} E[ x_k^2 | b_{l-1} < x_k \leq b_l ] \Bigg) \nonumber\\ = \; & \Phi \left( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus \{ k, k+1 \}} \beta_j x_j + \beta_k \bar{x}_{km} + \beta_{k+1} \overline{x^2}_{km} \right) \label{eq:probQuadEff} \\ & - \Phi \left( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus \{ k, k+1 \}} \beta_j x_j + \beta_k \bar{x}_{kl} + \beta_{k+1} \overline{x^2}_{kl} \right), \nonumber \end{align} where $\bar{x}_{km} \equiv E[ x_k | b_{m-1} < x_k \leq b_m ]$ and $\overline{x^2}_{km} \equiv E[ x_k^2 | b_{m-1} < x_k \leq b_m ]$ $\; \forall \; m = 1, \ldots, M$ remain the same as outlined in section~\ref{sec:linEffInt}.% \footnote{ See footnote~\ref{foot:PhiNonLin} regarding the approximate derivation from equation~(\ref{eq:probQuadEffStart}) to equation~(\ref{eq:probQuadEffMiddle}). } If the values of $\bar{x}_{km}$ and $\overline{x^2}_{km}$ are unknown, they can again be approximated by equations~(\ref{eq:linEffectXBar}) and~(\ref{eq:linEffectXSquaredBar}), respectively. In order to calculate the approximate standard error of $E_{k, ml}$, we again apply the Delta method: \begin{equation} \se( E_{k, ml}) = \sqrt{ \left( \frac{ \partial E_{k, ml}} { \partial \mathbf{\beta}} \right)^{\top} \cdot \var( \mathbf{\beta}) \cdot \frac{ \partial E_{k, ml}} { \partial \mathbf{\beta }}}, \label{eq:probQuadEffSE} \end{equation} where the elements of the gradient vector $\partial E_{k, ml} / \partial \mathbf{\beta}$ are: \begin{align} \frac{\partial E_{k, ml}}{\partial \beta_0} & = \phi_m - \phi_l \label{eq:intEffProb} \\ \frac{\partial E_{k, ml}}{\partial \beta_j} & = \left( \phi_m - \phi_l \right) x_j \; \forall \; j \in \{ 1, \ldots , k-1, k+2, \ldots , K \} \\ \frac{\partial E_{k, ml}}{\partial \beta_k} & = \phi_m \; \bar{x}_{km} - \phi_l \; \bar{x}_{kl} \\ \frac{\partial E_{k, ml}}{\partial \beta_{k+1}} & = \phi_m \; \overline{x^2}_{km} - \phi_l \; \overline{x^2}_{kl} \end{align} with $\phi_m \equiv \phi\left( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus \{ k , k+1 \} } \beta_j x_j + \beta_k \bar{x}_{km} + \beta_{k+1} \overline{x^2}_{km} \right) \; \forall m = \; 1, \ldots, M$. If the covariances between the estimated parameters are unknown, we can approximate equation~(\ref{eq:probQuadEffSE}) by assuming that the off-diagonal elements of the variance covariance matrix are zero. Using helper functions \code{EXSquared} (equation~\ref{eq:linEffectXSquaredBar}), \code{elaIntBounds} (checking arguments \code{refBound} and \code{intBound}), \code{checkXPos} (checking argument \code{xPos}, see appendix~\ref{sec:helperFunctions}), and \code{checkXBeta} (checking $\beta_0 + \sum_{j=1}^K \beta_j x_j$ for plausible values, see appendix~\ref{sec:helperFunctions}), the following function calculates the effect and its standard error according to equations~(\ref{eq:probEff}), (\ref{eq:probQuadEff}), and (\ref{eq:probQuadEffSE}): << >>= probEffInt <- function( allCoef, allXVal, xPos, refBound, intBound, allCoefSE = rep( NA, length( allCoef ) ) ){ # number of coefficients nCoef <- length( allCoef ) # check arguments if( length( allXVal ) != nCoef ){ stop( "argument 'allCoef' and 'allXVal' must have the same length" ) } if( length( allCoefSE ) != nCoef ){ stop( "argument 'allCoef' and 'allCoefSE' must have the same length" ) } checkXPos( xPos, minLength = 1, maxLength = 2, minVal = 1, maxVal = nCoef ) refBound <- elaIntBounds( refBound, 1, argName = "refBound" ) intBound <- elaIntBounds( intBound, 1, argName = "intBound" ) if( any( !is.na( allXVal[ xPos ] ) ) ) { allXVal[ xPos ] <- NA warning( "values of argument 'allXVal[ xPos ]' are ignored", " (set these values to 'NA' to avoid this warning)" ) } # calculate xBars intX <- mean( intBound ) refX <- mean( refBound ) if( length( xPos ) == 2 ) { intX <- c( intX, EXSquared( intBound[1], intBound[2] ) ) refX <- c( refX, EXSquared( refBound[1], refBound[2] ) ) } if( length( intX ) != length( xPos ) || length( refX ) != length( xPos ) ) { stop( "internal error: 'intX' or 'refX' does not have the expected length" ) } # define X' * beta intXbeta <- sum( allCoef * replace( allXVal, xPos, intX ) ) refXbeta <- sum( allCoef * replace( allXVal, xPos, refX ) ) checkXBeta( c( intXbeta, refXbeta ) ) # effect E_{k,ml} eff <- pnorm( intXbeta ) - pnorm( refXbeta ) # partial derivative of E_{k,ml} w.r.t. all estimated coefficients derivCoef <- rep( NA, nCoef ) derivCoef[ -xPos ] = ( dnorm( intXbeta ) - dnorm( refXbeta ) ) * allXVal[ -xPos ] derivCoef[ xPos ] = dnorm( intXbeta ) * intX - dnorm( refXbeta ) * refX # variance covariance of the coefficients (covariances set to zero) vcovCoef <- diag( allCoefSE^2 ) # approximate standard error of the effect effSE <- drop( sqrt( t( derivCoef ) %*% vcovCoef %*% derivCoef ) ) # object to be returned result <- c( effect = eff, stdEr = effSE ) return( result ) } @ where argument \code{allCoef} $= ( \beta_0, \ldots, \beta_K )^\top$, or alternatively for an ordered probit model \code{allCoef} $= ( -\mu_{p^* - 1}, \beta_1, \ldots, \beta_K )^\top$, is a vector containing all coefficients; argument \code{allXVal} $= ( 1, x_1, \ldots, x_K )^\top$ is a vector containing all values of the corresponding explanatory variables; argument \code{xPos} = $k+1$ or $( k+1 , k+2 )^{\top}$ is a scalar or vector with two elements that indicates the position of the variable of interest and its coefficient and eventually the position of the squared variable of interest and its coefficient in vectors \code{allXVal} and \code{allCoef}; argument \code{refBound} $= ( b_{l-1}, b_l )$ indicates the boundaries of the reference interval; argument \code{intBound} $= ( b_{m-1}, b_m )$ indicates the boundaries of the interval of interest; and optional argument \code{allCoefSE} $= ( \se(\beta_0), \ldots, \se(\beta_K) )^\top$, or alternatively for ordered probit models \code{allCoefSE} $= ( \se(\mu_{p^* - 1}), \se(\beta_1), \ldots, \se(\beta_K) )^\top$, is a vector of standard errors. Please note that the value of $x_k$ and% ---in case of an additional quadratic term---% also the value of $x_{k+1} = x_k^2$ in argument \code{allXVal}, i.e., \code{allXVal[xPos]}, are not used in the calculations of the effect~$E_{k, ml}$ or in the calculation of its standard error and, thus, the value of $x_k$ and---if present---also the value of $x_{k+1} = x_k^2$ should be set to~\code{NA}. << >>= # Example eff4a <- probEffInt( allCoef = c( 0.33, 0.22, 0.05, 0.6 ), allXVal = c( 1, NA, 0.16, 0.13 ), xPos = 2, refBound = c( 8, 12 ), intBound = c( 13, 15 )) eff4a eff4b <- probEffInt( allCoef = c( 0.33, 0.22, 0.05, 0.6 ), allXVal = c( 1, NA, 0.16, 0.13 ), xPos = 2, refBound = c( 8, 12 ), intBound = c( 13, 15 ), allCoefSE = c( 0.003, 0.045, 0.007, 0.009 ) ) eff4b # Example eff5a <- probEffInt( allCoef = c( 0.33, 0.22, 0.05, 0.006 ), allXVal = c( 1, NA, NA, 0.13 ), xPos = c( 2, 3 ), refBound = c( 8, 12 ), intBound = c( 13, 15 )) eff5a eff5b <- probEffInt( allCoef = c( 0.33, 0.22, 0.05, 0.006 ), allXVal = c( 1, NA, NA, 0.13 ), xPos = c( 2, 3 ), refBound = c( 8, 12 ), intBound = c( 13, 15 ), allCoefSE = c( 0.003, 0.045, 0.007, 0.009 ) ) eff5b @ %========================================================================================== \subsection{Grouping and re-basing categorical variables} \label{sec:2.4} We consider a regression model \begin{align} \Pr( y = 1 | x ) &= \Phi\left( \beta_0 + \sum_{j\in\{1, \ldots, K \} \setminus k} \beta_j x_j + \sum_{m \in \{1,...,M \} \setminus m^*} \delta_m D_m \right) \label{eq:probInterv}\\ D_m &= \begin{cases} 1 & \text{ if } x_k \in c_m \\ 0 & \text{ otherwise} \end{cases} \quad \forall \; m = 1, \ldots , M. \end{align} Like in section~\ref{sec:lingroup}, the $k$th explanatory variable is coded into $M$ mutually exclusive categories $c_1, \ldots, c_M$, with $c_m \cap c_l = \emptyset \; \forall \; m \neq l$, and $D_1, \ldots, D_M$ the corresponding dummy variables. In the case of a probit regression, equation~(\ref{eq:linEffGr}) modifies to: \begin{align} E_{k,lr} = \; & E[ y | x_k \in c_l^* ] - E[ y | x_k \in c_r^* ] \\ = \; & \Phi \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m E[ D_m | x_k \in c_l^* ] \right) \\ & - \Phi \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^{M} \delta_m E[ D_m | x_k \in c_r^* ] \right) \nonumber \\ = \; & \Phi \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m D_{ml} \right) \label{eq:probEffGr} \\ & - \Phi \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^{M} \delta_m D_{mr} \right) \nonumber \end{align} where $D_{mn}$ is defined as in equation~(\ref{eq:linEffGrD}). In the case of an ordered probit regression, $\beta_0$ in equation~(\ref{eq:probEffGr}) is again replaced by $-\mu_{p^* - 1}$. In order to calculate the approximate standard error of the new effect $E_{k,lr}$, we again apply the Delta method: \begin{equation} \se( E_{k,lr} ) = \sqrt{ \left( \frac{ \partial E_{k,lr}} { \partial \begin{pmatrix} \beta \\ \delta \end{pmatrix}} \right)^\top \cdot \var \begin{pmatrix} \beta \\ \delta \end{pmatrix} \cdot \frac{ \partial E_{k,lr}} { \partial \begin{pmatrix} \beta \\ \delta \end{pmatrix}}} , \label{eq:24Delta} \end{equation} with the partial derivatives defined as: \begin{align} \frac{\partial E_{k,lr}}{\partial \beta_0} &= \phi_{ml} (\cdot) - \phi_{mr} (\cdot) \\ \frac{\partial E_{k,lr}}{\partial \beta_j} &= ( \phi_{ml} (\cdot) - \phi_{mr} (\cdot)) \cdot \bar{x}_j \; \forall \; j \in \{ 1, \ldots , K \} \setminus k \\ \frac{\partial E_{k,lr}}{\partial \delta_m} &= \phi_{ml} (\cdot) D_{ml} - \phi_{mr} (\cdot) D_{mr} \; \forall \; m = 1, \ldots, M. \end{align} with \begin{align} \phi_{ml} (\cdot) & \equiv \phi \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m D_{ml} \right) \\ \phi_{mr} (\cdot) & \equiv \phi \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m D_{mr} \right) \label{eq:24grad} \end{align} Using the same example as in section~\ref{sec:lingroup}, equation~(\ref{eq:effectRB}) modifies to \begin{align} E_{k,\{1,2\}\{3,4\}} = \; & \Phi\left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \frac{ \delta_1 \cdot s_1 }{ s_1 + s_2 } + \frac{ \delta_2 \cdot s_2 }{ s_1 + s_2 } + \frac{ \delta_3 \cdot 0 }{ s_1 + s_2 } + \frac{ \delta_4 \cdot 0 }{ s_1 + s_2 } + \frac{ \delta_5 \cdot 0 }{ s_1 + s_2 } \right) \\ & - \Phi\left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \frac{ \delta_1 \cdot 0 }{ s_3 + s_4 } + \frac{ \delta_2 \cdot 0 }{ s_3 + s_4 } + \frac{ \delta_3 \cdot s_3 }{ s_3 + s_4 } + \frac{ \delta_4 \cdot s_4 }{ s_3 + s_4 } + \frac{ \delta_5 \cdot 0 }{ s_3 + s_4 } \right) \nonumber \\ = \; & \Phi\left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \frac{ \delta_1 \cdot s_1 + \delta_2 \cdot s_2 }{ s_1 + s_2 } \right) \\ & - \Phi\left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \frac{ \delta_3 \cdot s_3 + \delta_4 \cdot s_4 }{ s_3 + s_4 } \right) \nonumber \end{align} The following function calculates the effect $E_{k,lr}$ and its standard error in accordance with equations~(\ref{eq:probEffGr}) and (\ref{eq:24Delta}) to (\ref{eq:24grad}): << >>= probEffGr <- function( allCoef, allXVal, xPos, Group, allCoefSE = rep( NA, length( allCoef ) ) ){ nCoef <- length( allCoef ) xShares <- allXVal[ xPos ] xCoef <- allCoef[ xPos ] if( sum( xShares ) > 1 ){ stop( "the shares in argument 'xShares' sum up to a value larger than 1" ) } if( length( xCoef ) != length( xShares ) ){ stop( "arguments 'xCoef' and 'xShares' must have the same length" ) } if( length( xCoef ) != length( Group ) ){ stop( "arguments 'xCoef' and 'Group' must have the same length" ) } if( ! all( Group %in% c( -1, 0, 1 ) ) ){ stop( "all elements of argument 'Group' must be -1, 0, or 1" ) } # D_mr DRef <- sum( xCoef[ Group == -1 ] * xShares[ Group == -1 ]) / sum( xShares[ Group == -1 ] ) XBetaRef <- sum( allCoef[ -xPos ] * allXVal[ -xPos ]) + DRef # D_ml DEffect <- sum( xCoef[ Group == 1 ] * xShares[ Group == 1 ]) / sum( xShares[ Group == 1 ] ) XBetaEffect <- sum( allCoef[ -xPos ] * allXVal[ -xPos ]) + DEffect # effect effeG <- pnorm( XBetaEffect ) - pnorm( XBetaRef ) # partial derivative of E_{k,ml} w.r.t. all estimated coefficients derivCoef <- rep( NA, nCoef ) derivCoef[ -xPos ] = ( dnorm( XBetaEffect ) - dnorm( XBetaRef ) ) * allXVal[ -xPos ] derivCoef[ xPos ] = dnorm( XBetaEffect ) * DEffect - dnorm( XBetaRef ) * DRef # variance covariance of the coefficients (covariances set to zero) vcovCoef <- diag( allCoefSE^2 ) # approximate standard error of the effect effeGSE <- drop( sqrt( t( derivCoef ) %*% vcovCoef %*% derivCoef ) ) result <- c( effect = effeG, stdEr = effeGSE ) return( result ) } @ where argument \code{allCoef} $= (\beta_0, \ldots, \beta_{k-1}, \delta_1, \ldots, \delta_M, \beta_{k+1}, \ldots, \beta_K )^\top$ or alternatively \code{allCoef} $= ( -\mu_{p^* -1}, \beta_1, \ldots, \beta_{k-1}, \delta_1, \ldots, \delta_M, \beta_{k+1}, \ldots, \beta_K )^\top$ are the coefficients of the probit regression or ordered probit regression, respectively, \textbf{with the coefficient of the reference group of the $k$th variable set to zero}, i.e., $\delta_{m^*} \equiv 0$; argument \code{allXVal} $= (x_1, \ldots, x_{k-1}, s_1, \ldots, s_M, x_{k+1}, \ldots, x_K )^\top$ are the mean values of the respective explanatory variables and the shares of each of the categories of the $k$th variable, \textbf{which includes the share of the reference group of the $k$th variable} ($s_{m^*}$); argument \code{xPos} is a vector of $M$~elements indicating the positions of the coefficients~$\delta_1, \ldots, \delta_M$ in argument \code{allCoef} and, equally, the positions of the shares~$s_1, \ldots, s_M$ in argument \code{allXVal}; argument \code{Group} $= ( D_{1l} - D_{1r}, \ldots , D_{Ml} - D_{Mr} )^{\top}$ is a vector of $M$~elements consisting of $-1$s, $0$s, and $1$s, where a $-1$ indicates that the category belongs to the new reference category, a $1$ indicates that the category belongs to the new category of interest, and a zero indicates that the category belongs to neither; finally optional argument \code{allCoefSE} $=( \se( \beta_0), \ldots, \se( \beta_{k-1}),$ $\se(\delta_1), \ldots, \se(\delta_M ),$ $\se( \beta_{k+1}), \ldots, \se( \beta_K) )^\top$, or alternatively in the case of an ordered probit model \code{allCoefSE} $=( \se( \mu_{m^* - 1} ), \se( \beta_1 ), \ldots, \se( \beta_{k-1}),$ $\se(\delta_1), \ldots, \se(\delta_M ),$ $\se( \beta_{k+1}), \ldots, \se( \beta_K) )^\top$ are the standard errors of all coefficients of the probit regression, \textbf{where the standard error for the reference group is included as zero}, i.e., $\se(\delta_{m^*}) \equiv 0$. Elements of arguments \code{allCoef}, \code{allXVal}, \code{xPos}, \code{Group}, and \code{allCoefSE} that belong to a category that is neither in the reference category nor in the category of interest, i.e., categories~$m$ with $D_{mr} = D_{ml} = 0$, can be omitted but the ommission must be consistent for all five arguments. << >>= # Example Eff6a <- probEffGr( allCoef = c( 0.28, 0.003, 0.0075, 0, -0.034, -0.005, 0.89, -1.2 ), allXVal = c( 1, 0.1, 0.3, 0.25, 0.15, 0.2, 2.34, 10.8 ), xPos = c( 2:6 ), Group = c( 0, -1, -1, 1, 1 ) ) Eff6a # Example Eff6b <- probEffGr( allCoef = c( 0.28, 0.003, 0.0075, 0, -0.034, -0.005, 0.89, -1.2 ), allXVal = c( 1, 0.1, 0.3, 0.25, 0.15, 0.2, 2.34, 10.8 ), xPos = c( 2:6 ), Group = c( 0, -1, -1, 1, 1 ), allCoefSE = c( 0.03, 0.0001, 0.005, 0, 0.01, 0.004, 0.05, 0.8 )) Eff6b # the same examples but without categories that are neither # in the new reference category nor in the new category of interest all.equal( Eff6a, probEffGr( allCoef = c( 0.28, 0.0075, 0, -0.034, -0.005, 0.89, -1.2 ), allXVal = c( 1, 0.3, 0.25, 0.15, 0.2, 2.34, 10.8 ), xPos = c( 2:5 ), Group = c( -1, -1, 1, 1 ) ) ) all.equal( Eff6b, probEffGr( allCoef = c( 0.28, 0.0075, 0, -0.034, -0.005, 0.89, -1.2 ), allXVal = c( 1, 0.3, 0.25, 0.15, 0.2, 2.34, 10.8 ), xPos = c( 2:5 ), Group = c( -1, -1, 1, 1 ), allCoefSE = c( 0.03, 0.005, 0, 0.01, 0.004, 0.05, 0.8 ) ) ) @ %========================================================================================== \section{Binary logit model, multinomial logit model, conditional logit model, and nested logit model} \subsection{Semi-elasticities from continuous explanatory variables (linear and quadratic)} \label{seq:logit} Binary logit model with continuous explanatory variables: \begin{equation} \Pr( y = 1 | x ) = \frac{ \exp{( \beta_0 + \sum_{k=1}^K \beta_k x_k )}} { 1 + \exp{( \beta_0 + \sum_{k=1}^K \beta_k x_k )}}, \label{eq:logit} \end{equation} where $y \in \{0, 1 \}$ is a binary dependent variable, $x = (x_1, \ldots, x_K )^{\top}$ is a vector of $K$ continuous explanatory variables, and $\beta = ( \beta_0, \ldots, \beta_K )^{\top}$ is a vector of $K + 1$ unknown coefficients. The semi-elasticity of the $k$th (continuous) explanatory variable for the binary logit model is: \begin{equation} \epsilon_k = \frac{\partial \Pr( y = 1 )}{\partial \bar{x}_k } \cdot \bar{x}_k = \frac{ \exp{( \beta_0 + \sum_{k=1}^K \beta_k \bar{x}_k )}} {( 1 + \exp{( \beta_0 + \sum_{k=1}^K \beta_k \bar{x}_k )})^2 } \cdot \beta_k \bar{x}_k \label{eq:logEla1} \end{equation} The interpretation is equal to the one described in section~\ref{sec:LPMcont}. If $x_k$ also enters equation~(\ref{eq:logit}) in quadratic form the semi-elasticity modifies to \begin{equation} \epsilon_k = \frac{ \exp{( \beta_0 + \sum_{k \in \{1,\ldots,K \}\setminus k+1 } \beta_k \bar{x}_k + \beta_{k+1} \bar{x}_k^2 )}}{ \left( 1 + \exp{( \beta_0 + \sum_{k \in \{1,\ldots, K\}\setminus k+1 } \beta_k \bar{x}_k + \beta_{k+1} \bar{x}_k^2 )} \right)^2} \cdot ( \beta_k \bar{x}_k + 2 \beta_{k+1} \bar{x}_k^2 ) \label{eq:logEla2} \end{equation} To calculate the (approximate) standard error of $\epsilon_k$, we again use the simplified Delta method as described in section~\ref{sec:probEla}, i.e., we assume the covariances to be zero and $\frac{ \exp{( \beta_0 + \sum_{k=1}^K \beta_k \bar{x}_k )}} {( 1 + \exp{( \beta_0 + \sum_{k=1}^K \beta_k \bar{x}_k )})^2 }$ to be a constant. The standard error then calculates as \begin{align} \se(\epsilon_k) &= \sqrt{ \frac{ \exp{( \beta_0 + \sum_{k=1}^K \beta_k \bar{x}_k )}} {( 1 + \exp{( \beta_0 + \sum_{k=1}^K \beta_k \bar{x}_k )})^2 } \cdot \bar{x}_k \cdot \var( \beta_k ) \cdot \frac{ \exp{( \beta_0 + \sum_{k=1}^K \beta_k \bar{x}_k )}} {( 1 + \exp{( \beta_0 + \sum_{k=1}^K \beta_k \bar{x}_k )})^2 } \cdot \bar{x}_k } \\ &= \frac{ \exp{( \beta_0 + \sum_{k=1}^K \beta_k \bar{x}_k )}} {( 1 + \exp{( \beta_0 + \sum_{k=1}^K \beta_k \bar{x}_k )})^2 } \cdot \bar{x}_k \cdot \se( \beta_k ) \label{eq:logEla1SE} \end{align} and in case of a quadratic covariate $\bar{x}_k^2$ \begin{align} \frac{\partial \epsilon_k}{\partial \begin{pmatrix} \beta_k \\ \beta_{k+1} \end{pmatrix}} &= \begin{pmatrix} \dfrac{ \exp{( \beta_0 + \sum_{k \in \{1,\ldots,K \}\setminus k+1 } \beta_k \bar{x}_k + \beta_{k+1} \bar{x}_k^2 )}}{ \left( 1 + \exp{( \beta_0 + \sum_{k \in \{1,\ldots, K\}\setminus k+1 } \beta_k \bar{x}_k + \beta_{k+1} \bar{x}_k^2 )} \right)^2} \cdot \bar{x}_k \\ \\ \dfrac{ \exp{( \beta_0 + \sum_{k \in \{1,\ldots,K \}\setminus k+1 } \beta_k \bar{x}_k + \beta_{k+1} \bar{x}_k^2 )}}{ \left( 1 + \exp{( \beta_0 + \sum_{k \in \{1,\ldots, K\}\setminus k+1 } \beta_k \bar{x}_k + \beta_{k+1} \bar{x}_k^2 )} \right)^2} \cdot 2\bar{x}_k^2 \end{pmatrix} \\ \se( \epsilon_k) &= \sqrt{ \left( \frac{\partial \epsilon_k}{\partial \begin{pmatrix} \beta_k \\ \beta_{k+1} \end{pmatrix}}\right)^{\top} \cdot \var \begin{pmatrix} \beta_k \\ \beta_{k+1} \end{pmatrix} \cdot \frac{\partial \epsilon_k}{\partial \begin{pmatrix} \beta_k \\ \beta_{k+1} \end{pmatrix}}} \label{eq:logEla2SE} \end{align} %================================================================================================================ The multinomial logit model with continuous explanatory variables: \begin{equation} \Pr( y = p | x ) = \pi_p = \frac{ \exp( \beta_{0,p} + \sum_{k=1}^K \beta_{k,p} x_k ) )} { 1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp( \beta_{0,p} + \sum_{k=1}^K \beta_{k,p} x_k ) )} \label{eq:MNLogit} \end{equation} where $y \in \{1, \ldots, P \}$ is a categorical dependent variable of which the base category, $p*$, is set to zero, $x = (x_1, \ldots, x_K )^{\top}$ is a vector of $K$ continuous explanatory variables describing the characteristics of the choice taker and $\beta = ( \beta_{0, p}, \ldots, \beta_{K, p} )^{\top}$ is a vector of $K + 1$ unknown coefficients. The semi-elasticity of the $k$th (continuous) explanatory variable for the multinomial logit model is: \begin{align} \epsilon_{k,p} &= \frac{ \partial \Pr( y = p )}{ \partial \bar{x}_k} \cdot \bar{x}_k \nonumber \\ &= \pi_p \left( \frac{ \beta_{k,p}}{[ 1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp( \beta_{0,p} + \sum_{k=1}^K \beta_{k,p} x_k )]} + \sum_{o \in \{1, \ldots, P \} \setminus \{ p, p^* \}} ( \beta_{k,p} - \beta_{k,o}) \pi_o \right) \cdot \bar{x}_k \label{eq:MNLogitEla} \end{align} This semi-elasticity can be interpreted as: if the $k,p$th explanatory variable inceases by one percent, the probability that $y$ belongs to category $p$ instead of the base category, increases by $\epsilon_{k,p}$ percentage points. If $x_k$ also enters equation~(\ref{eq:MNLogit}) in quadratic form, \begin{equation} \Pr( y = p | x ) = \pi_p = \frac{ \exp( \beta_{0,p} + \sum_{k \in \{1,\ldots,K \}\setminus k+1 } \beta_{k,p} \bar{x}_k + \beta_{k+1,p} \bar{x}_k^2 ) )} { 1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp( \beta_{0,p} + \sum_{k \in \{1,\ldots,K \}\setminus k+1 } \beta_{k,p} \bar{x}_k + \beta_{k+1,p} \bar{x}_k^2 ) )} \label{eq:MNLogitQuad} \end{equation} equation~(\ref{eq:MNLogitEla}) modifies to \begin{align} \epsilon_{k,p} =& \frac{ \partial \Pr( y = p )}{ \partial \bar{x}_k} \cdot \bar{x}_k \nonumber \\ =& \pi_p \left( \frac{ ( \beta_{k,p} + 2 \beta_{k+1,p} \bar{x}_k )} {[ 1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp( \beta_{0,p} + \sum_{k \in \{1,\ldots,K \}\setminus k+1 } \beta_{k,p} \bar{x}_k + \beta_{k+1,p} \bar{x}_k^2)]} \right. \\ & \left. + \sum_{o \in \{1, \ldots, P \} \setminus \{ p, p^* \}} (( \beta_{k,p} + 2 \beta_{k+1,p} \bar{x}_k ) - ( \beta_{k,o} + 2 \beta_{k+1,0} \bar{x}_k )) \pi_o \right) \cdot \bar{x}_k \nonumber \label{eq:MNLogitEla2} \end{align} To calculate the (approximate) standard error $\epsilon_{k,p}$, we use a simplified Delta method, where we again set all covariances between the $\beta$s to zero and assume for the most part fixed factors for the remaining terms. Hence, \begin{align} \frac{\partial \epsilon_{k,p}}{\partial \beta_{k,p}} =& \frac{ \exp( \beta_{0,p} + \sum_{k=1}^K \beta_{k,p} x_k ) )} { [ 1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp( \beta_{0,p} + \sum_{k=1}^K \beta_{k,p} x_k ) ]^2 } \bar{x}_k \\ & + \pi_p \sum_{o \in \{1, \ldots, P \} \setminus \{ p, p^* \}} ( \beta_{k,p} - \beta_{k,o} ) \pi_o \cdot \bar{x}_k \nonumber \end{align} \begin{align} \se( \epsilon_{k,p}) &= \sqrt{ \left( \frac{\partial \epsilon_{k,p}}{\partial \beta_{k,p}} \right)^\top \var( \beta_{k,p}) \frac{\partial \epsilon_{k,p}}{\partial \beta_{k,p}} } \nonumber \\ &= \sqrt{ \frac{\partial \epsilon_{k,p}}{\partial \beta_{k,p}}} \se( \beta_{k,p}) \label{eq:eq:MNlogElaSE} \end{align} and in case of a quadratic covariate $\bar{x}_k^2$ \begin{align} \frac{\partial \epsilon_{k,p}}{\partial \begin{pmatrix} \beta_{k,p} \\ \beta_{k+1,p} \end{pmatrix}} &= \begin{pmatrix} \dfrac{ \exp( \beta_{0,p} + \sum_{k \in \{1,\ldots,K \}\setminus k+1 } \beta_{k,p} \bar{x}_k + \beta_{k+1,p} \bar{x}_k^2 ) ) )} { [ 1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp( \beta_{0,p} + \sum_{k \in \{1,\ldots,K \}\setminus k+1 } \beta_{k,p} \bar{x}_k + \beta_{k+1,p} \bar{x}_k^2 ) ]^2 } \bar{x}_k \\ \\ + \pi_p \sum_{o \in \{1, \ldots, P \} \setminus \{ p, p^* \}} ( ( \beta_{k,p} + 2 \beta_{k+1,p} \bar{x}_k ) - ( \beta_{k,o} + 2 \beta_{k+1,0} \bar{x}_k ) ) \pi_o \cdot \bar{x}_k \\ \\ \\ \dfrac{ \exp( \beta_{0,p} + \sum_{k \in \{1,\ldots,K \}\setminus k+1 } \beta_{k,p} \bar{x}_k + \beta_{k+1,p} \bar{x}_k^2 ) ) )} { [ 1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp( \beta_{0,p} + \sum_{k \in \{1,\ldots,K \}\setminus k+1 } \beta_{k,p} \bar{x}_k + \beta_{k+1,p} \bar{x}_k^2 ) ]^2 } 2 \bar{x}_k^2 \\ \\ + \pi_p \sum_{o \in \{1, \ldots, P \} \setminus \{ p, p^* \}} ( ( \beta_{k,p} + 2 \beta_{k+1,p} \bar{x}_k ) - ( \beta_{k,o} + 2 \beta_{k+1,0} \bar{x}_k ) ) \pi_o \cdot 2 \bar{x}_k^2 \end{pmatrix} \\ \se( \epsilon_{k,p}) &= \sqrt{ \left( \frac{\partial \epsilon_{k,p}} {\partial \begin{pmatrix} \beta_{k,p} \\ \beta_{k+1,p} \end{pmatrix}} \right)^\top \cdot \var \begin{pmatrix} \beta_{k,p} \\ \beta_{k+1,p} \end{pmatrix} \cdot \frac{\partial \epsilon_{k,p}} {\partial \begin{pmatrix} \beta_{k,p} \\ \beta_{k+1,p} \end{pmatrix}}} \label{eq:MNlogEla2SE} \end{align} %============================================================================================================ Conditional logit model with continuous explanatory variables: \begin{equation} \Pr( y = p | z ) = \pi_p = \frac{ \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,p} ) ) } { \sum_{p \in C } \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,p} ) ) } \label{eq:MNLogit} \end{equation} where $y \in \{1, \ldots, P \}$ is a categorical dependent variable, $z = (z_{1,p}, \ldots, z_{T,p})^{\top}$ is a vector of $T$ continuous explanatory variables describing the characteristics of each of the $P$ alternatives in the choice set $C$, and $\gamma = ( \gamma_0, \ldots, \gamma_T )^{\top}$ is a vector of $T + 1$ unknown coefficients. The semi-elasticity of the $t$th (continuous) explanatory variable for the conditional logit model is: \begin{align} \epsilon_{t,p} &= \frac{\partial \Pr( y = p )}{\partial \bar{z}_{t,p} } \cdot \bar{z}_{t,p} \nonumber \\ &= ( \pi_p - \pi_p^2 ) \cdot \gamma_t \bar{z}_{t,p} \label{eq:elaCondLog} \end{align} This semi-elasticity can be interpreted as: if the $t$th explanatory variable inceases by one percent, the probability that $y$ belongs to category $j$ increases by $\epsilon_{t,p}$ percentage points. If $z_{t,p}$ is also included in quadratic form, equation~(\ref{eq:elaCondLog}) modifies to \begin{align} \epsilon_{t,p} &= \frac{\partial \Pr( y = p )}{\partial \bar{z}_{t,p} } \cdot \bar{z}_{t,p} \nonumber \\ &= ( \pi_p - \pi_p^2 ) \cdot ( \gamma_t \bar{z}_{t,p} + 2 \gamma_{t+1} \bar{z}_{t,p}^2 ) \label{eq:elaCondLog2} \end{align} In order to calculate approximate standard errors in the case of the conditional logit regression, equation~(\ref{eq:logEla1SE}) modifies to \begin{equation} \se( \epsilon_{t,p}) = ( \pi_p - \pi_p^2 ) z_{t,p} \cdot \se( \gamma_t ) %\se( \epsilon_{t,p}) = (( 1 - 3 \pi_p + 2 \pi_p^2 ) \gamma_t z_{t,p} - \pi_p + 1 ) \pi_p z_{t,p} \cdot \se( \gamma_t ) \label{eq:CondLogElaSE} \end{equation} and in case of a quadratic covariate $\bar{z}_t^2$ \begin{align} \frac{\partial \epsilon_{t,p}}{\partial \begin{pmatrix} \gamma_{t} \\ \gamma_{t+1} \end{pmatrix}} &= \begin{pmatrix} ( \pi_p - \pi_p^2 ) z_{t,p} \\ \\ ( \pi_p - \pi_p^2 ) 2 z_{t,p}^2 \end{pmatrix} \\ \se( \epsilon_{t,p}) &= \sqrt{ \left( \frac{\partial \epsilon_{t,p}} {\partial \begin{pmatrix} \gamma_{t} \\ \gamma_{t+1} \end{pmatrix}} \right)^\top \cdot \var\begin{pmatrix} \gamma_{t} \\ \gamma_{t+1} \end{pmatrix} \cdot \frac{\partial \epsilon_{t,p}} {\partial \begin{pmatrix} \gamma_{t} \\ \gamma_{t+1} \end{pmatrix}}} \label{eq:CondLogEla2SE} \end{align} %======================================================================================================= For the nested logit model with continuous explanatory variables at both the branch and the twig level,% \footnote{We restrict our analysis on a nesting structure with two levels, branches and twigs, as these seem most common in the applied empirical work. For nested logits with higher levelled nesting structures, the formulas must be modified accordingly. However, the reader should not that an analytical solution for nested logits with higher levelled nesting structures might become very complicated and computational solutions might be preferred.} we have to differentiate between the probability of choice of the $o$th branch, $B_o$, \begin{align} \Pr( y \in B_o | x ) = \pi_o =& \; \frac{ \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_o ) } { \sum_{o = 1}^O \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_o )} \label{eq:semelaBranch} \\ \text{ with } \;\;\; IV_o =& \; \ln \sum_{p \in B_o} \exp \left( \beta_0/\lambda_o + \sum_{k=1}^K \beta_k/\lambda_o x_{k,p} \right) \label{eq:IVo} \end{align} and the combined probablity of choosing the $p$th alternative in the $o$th branch \begin{align} \Pr( y = p, y \in B_o ) =& \; \Pr( y=p | x, y \in B_o ) \cdot \Pr( y \in B_o | x ) = \pi_p \cdot \pi_o \nonumber \\[1em] =& \; \frac{ \exp( \beta_0/\lambda_o + \sum_{k=1}^K \beta_k/\lambda_o x_{k,p} )} { \sum_{p \in B_o} \exp( \beta_0/\lambda_o + \sum_{k=1}^K \beta_k/\lambda_o x_{k,p} )} \cdot \frac{ \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_o )} { \sum_{o = 1}^O \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_o )} \end{align} where $z_o = ( 1, z_{1,o}, \ldots, z_{T,o})^\top$ are the branch specific sets of characterics (which can be empty sets), $\gamma = (\gamma_0, \ldots, \gamma_T )^\top$ are the coefficients of the branch specific characteristics, $x_p = ( 1, x_{1,p}, \ldots, x_{K,p})^\top$ are variables describing the characteristics of the $P_o$ alternatives of $B_o$ and/or of the choice maker, and $\beta = ( \beta_0, \ldots, \beta_K )^\top$ are the corresponding coefficients, and $\lambda = ( \lambda_1, \ldots, \lambda_O)^\top$ are the parameters of the correlation coefficients ($1-\lambda_o$) between the $P_o$ branch specific alternatives. At the branch level, the semi-elasticity for a variable $z_{t,o}$ can be derived as: \begin{equation} \epsilon_{t,o} = ( \pi_o - \pi_o^2) \cdot \gamma_t \bar{z}_{t,o} \end{equation} and if $z_{t,o}$ enters the equation in quadratic form: \begin{equation} \epsilon_{t,o} = ( \pi_o - \pi_o^2) \cdot ( \gamma_t \bar{z}_{t,o} + 2 \gamma_{t+1} \bar{z}_{t,o}^2 ) \end{equation} whilst the semi-elasticity for $x_{k,p}$ wrt the combined probability is: \begin{equation} \epsilon_{k,p} = \left( \pi_o( \pi_p - \pi_p^2 ) \frac{1}{\lambda_o} + \pi_p^2 ( \pi_o - \pi_o^2 ) \lambda_o IV_o \right) \cdot \beta_k \bar{x}_{k,p} \end{equation} and for the quadratic form \begin{equation} \epsilon_{k,p} = \left( \pi_o( \pi_p - \pi_p^2 ) \frac{1}{\lambda_o} + \pi_p^2 ( \pi_o - \pi_o^2 ) \lambda_o IV_o \right) \cdot ( \beta_k \bar{x}_{k,p} + 2 \beta_{k+1} \bar{x}_{k,p}^2 ) \end{equation} In order to calculate approximate standard errors for the semi-elasticity at the branch level (equation~(\ref{eq:semelaBranch})), equation~(\ref{eq:logEla1SE}) modifies to \begin{equation} \se( \epsilon_{t,o}) = ( \pi_o - \pi_o^2 ) \cdot z_{t,o} \cdot \se( \gamma_t ) \label{eq:NestedLogElaSE} \end{equation} and in case of a quadratic covariate $\bar{z}_{t,o}^2$ \begin{align} \frac{\partial \epsilon_{t,o}}{\partial \begin{pmatrix} \gamma_{t} \\ \gamma_{t+1} \end{pmatrix}} &= \begin{pmatrix} ( \pi_o - \pi_o^2 ) z_{t,o} \\ \\ ( \pi_o - \pi_o^2 ) 2 z_{t,o}^2 \end{pmatrix} \\ \se( \epsilon_{t,o}) &= \sqrt{ \left( \frac{\partial \epsilon_{t,o}} {\partial \begin{pmatrix} \gamma_{t} \\ \gamma_{t+1} \end{pmatrix}} \right)^\top \cdot \var \begin{pmatrix} \gamma_{t} \\ \gamma_{t+1} \end{pmatrix} \cdot \frac{\partial \epsilon_{t,o}} {\partial \begin{pmatrix} \gamma_{t} \\ \gamma_{t+1} \end{pmatrix}}} \label{eq:NestedLogEla2SE} \end{align} And the standard error of $\epsilon_{k,p}$ can be calculated as \begin{equation} \se( \epsilon_{k,p}) = \left( \pi_o( \pi_p - \pi_p^2 ) \frac{1}{\lambda_o} + \pi_p^2 ( \pi_o - \pi_o^2 ) \lambda_o IV_o \right) \cdot x_{k,p} \cdot \se( \beta_k ) \label{eq:NestedLogElaSE2} \end{equation} and in case of a quadratic covariate $\bar{x}_{k,p}^2$ \begin{align} \frac{\partial \epsilon_{k,p}}{\partial \begin{pmatrix} \beta_{k} \\ \beta_{k+1} \end{pmatrix}} &= \begin{pmatrix} \left( \pi_o( \pi_p - \pi_p^2 ) \frac{1}{\lambda_o} + \pi_p^2 ( \pi_o - \pi_o^2 ) \lambda_o IV_o \right) x_{k,p} \\ \\ \left( \pi_o( \pi_p - \pi_p^2 ) \frac{1}{\lambda_o} + \pi_p^2 ( \pi_o - \pi_o^2 ) \lambda_o IV_o \right) 2 x_{k,p}^2 \end{pmatrix} \\ \se( \epsilon_{k,p}) &= \sqrt{ \left( \frac{\partial \epsilon_{k,p}} {\partial \begin{pmatrix} \beta_{k} \\ \beta_{k+1} \end{pmatrix}} \right)^\top \cdot \var \begin{pmatrix} \beta_{k} \\ \beta_{k+1} \end{pmatrix} \cdot \frac{\partial \epsilon_{k,p}} {\partial \begin{pmatrix} \beta_{k} \\ \beta_{k+1} \end{pmatrix}}} \label{eq:NestedLogEla2SE2} \end{align} Using the helper functions \code{checkXPos} (checking argument \code{xPos}, see appendix~\ref{sec:helperFunctions}) and \code{check-PKXBeta} (checking $\beta_0 + \sum_{j=1}^K \beta_j x_j$ for plausible values, see appendix~\ref{sec:helperFunctions}), the following function calculates the semi-elasticities and their respective standard errors as defined in equations~(\ref{eq:logEla1}), (\ref{eq:logEla2}), (\ref{eq:logEla1SE}), and (\ref{eq:logEla2SE}) for the binary logit function, the semi-elasticities and their respective approximate standard errors as defined in equations~(\ref{eq:MNLogitEla}), (\ref{eq:MNLogitEla2}), (\ref{eq:eq:MNlogElaSE}), and~(\ref{eq:MNlogEla2SE}) for the multi-nomial logit function, the semi-elasticities and their respective approcimate standard errors as defined in equations~(\ref{eq:elaCondLog}), (\ref{eq:elaCondLog2}), (\ref{eq:CondLogElaSE}), and~(\ref{eq:CondLogEla2SE}) for the conditional logit << >>= logEla <- function( allCoef, allCoefBra, allXVal, allXValBra, allCoefSE = rep( NA, length( allCoef ) ), xPos, yCat, yCatBra, lambda, method = "binary" ){ if( method == "binary" || method == "CondL" ){ nCoef <- length( allCoef ) # Checking standard errors if( nCoef != length( allCoefSE ) ) { stop( "arguments 'allCoef' and 'allCoefSE' must have the same length" ) } } else if( method == "MNL" ){ NCoef <- length( allCoef ) mCoef <- matrix( allCoef, nrow = length( allXVal ) ) nCoef <- dim( mCoef )[1] pCoef <- dim( mCoef )[2] # Checking standard errors if( NCoef != length( allCoefSE ) ) { stop( "arguments 'allCoef' and 'allCoefSE' must have the same length" ) } } else{ nCoef <- length( allCoef ) NCoef <- length( allCoefBra ) # Checking standard errors if( nCoef != length( allCoefSE ) ){ stop( "arguments 'allCoef' and 'allCoefSE' must have the same length" ) } } # Check position vector checkXPos( xPos, minLength = 1, maxLength = 2, minVal = 1, maxVal = nCoef ) # Check x values if( method == "binary" || method == "MNL" ){ if( nCoef != length( allXVal ) ) { stop( "arguments 'allCoef' and 'allXVal' must have the same length" ) } } else if( method == "CondL" ){ mXVal <- matrix( allXVal, nrow = nCoef ) nXVal <- dim( mXVal )[1] pXVal <- dim( mXVal )[2] if( nCoef != dim( mXVal )[1] ) { stop( "arguments 'allCoef' and 'allXVal' must have the same length" ) } } else{ mXValBra <- matrix( allXValBra, nrow = NCoef ) nXValBra <- dim( mXValBra )[1] pXValBra <- dim( mXValBra )[2] if( NCoef != nXValBra ) { stop( "arguments 'allCoefBra' and 'allXValBra' must have the same length" ) } O <- length( allXVal ) mXVal <- matrix( unlist( allXVal[ yCatBra ] ), nrow = nCoef ) nXVal <- dim( mXVal )[1] pXVal <- dim( mXVal )[2] if( nCoef != nXVal ) { stop( "arguments 'allCoef' and 'allXVal' must have the same length" ) } } # Identify coefficients of interest (kth/tth covariate) if( length( xPos ) == 2 ){ if( method == "binary" ){ xCoef <- allCoef[ xPos ] if( !isTRUE( all.equal( allXVal[xPos[2]], allXVal[xPos[1]]^2 ) ) ) { stop( "the value of 'allXVal[ xPos[2] ]' must be equal", "to the squared value of 'allXVal[ xPos[1] ]' " ) } } else if( method == "MNL" ){ xCoef <- mCoef[ xPos, ] if( !isTRUE( all.equal( allXVal[xPos[2]], allXVal[xPos[1]]^2 ) ) ) { stop( "the value of 'allXVal[ xPos[2] ]' must be equal", "to the squared value of 'allXVal[ xPos[1] ]' " ) } } else if( method == "CondL" ){ xCoef <- allCoef[ xPos ] for( p in 1:pXVal ){ if( !isTRUE( all.equal( mXVal[xPos[2], p], mXVal[xPos[1], p]^2 ) ) ) { stop( "the value of 'allXVal[ xPos[2] ]' must be equal", "to the squared value of 'allXVal[ xPos[1] ]' " ) } } } else{ xCoef <- allCoef[ xPos ] for( p in 1:pXVal ){ if( !isTRUE( all.equal( mXVal[xPos[2], p], mXVal[xPos[1], p]^2 ) ) ) { stop( "the value of 'allXVal[ xPos[2] ]' must be equal", "to the squared value of 'allXVal[ xPos[1] ]' " ) } } } } else if( length( xPos ) == 1 ) { if( method == "binary" || method == "CondL" ){ xCoef <- c( allCoef[ xPos ], 0 ) } else if( method == "MNL" ){ xCoef <- matrix( c( mCoef[ xPos, ], rep( 0, dim( mCoef )[ 2 ] ) ), nrow = 2, byrow = TRUE ) } else{ xCoef <- c( allCoef[ xPos ], 0 ) } } else { stop( "argument 'xPos' must be a scalar or a vector with two elements" ) } if( method == "binary" ){ xVal <- allXVal[ xPos[1] ] xBeta <- sum( allCoef * allXVal ) checkXBeta( xBeta ) dfun <- exp( xBeta )/( 1 + exp( xBeta ) )^2 semEla <- ( xCoef[ 1 ] + 2 * xCoef[ 2 ] * xVal ) * xVal * dfun } else if( method == "MNL" ){ #checkXBeta missing xVal <- allXVal[ xPos[1] ] xBeta <- colSums( mCoef * allXVal ) pfun <- rep( NA, length( xBeta )) term <- 0 for( i in 1:length( xBeta )){ pfun[i] <- exp( xBeta[i] )/( 1 + sum( exp( xBeta ) ) ) term <- term + ( ( xCoef[ 1, yCat ] + 2 * xCoef[ 2, yCat ] * xVal ) - ( xCoef[ 1, i ] + 2 * xCoef[ 2, i ] * xVal ) * pfun[i] ) } semEla <- xVal * pfun[ yCat ] * ( ( xCoef[ 1, yCat ] + 2 * xCoef[ 2, yCat ] * xVal )/ ( 1 + sum( exp( xBeta ))) + term ) dfun <- pfun[ yCat ] * ( 1/( 1 + sum( exp( xBeta ) ) ) + term ) } else if( method == "CondL" ){ #checkXBeta missing xVal <- rep( NA, pXVal ) for( p in 1:pXVal ){ xVal[p] <- mXVal[ xPos[ 1 ], p ] } xBeta <- colSums( allCoef * mXVal ) pfun <- exp( xBeta[ yCat ] )/( sum( exp( xBeta ) ) ) semEla <- ( xCoef[1] + 2 * xCoef[2] * xVal[ yCat ] ) * xVal[ yCat ] * ( pfun - pfun^2 ) } else{ #checkXBeta missing xVal <- rep( NA, pXVal ) for( p in 1:pXVal ){ xVal[p] <- mXVal[ xPos[ 1 ], p ] } coef <- matrix( NA, nrow = O, ncol = nCoef ) for( o in 1:O ){ coef[o, ] <- allCoef/lambda[o] } xBeta <- lapply( 1:O, function( i, m, v ){ colSums( m[[i]] * v[[i]] ) }, m=allXVal, v=coef ) IV <- unlist( lapply( 1:O, function( i, m ){ log( sum( exp( m[[i]] ) ) ) }, m=xBeta ) ) pfun <- exp( xBeta[[ yCatBra ]][ yCat ] )/ ( sum( exp( xBeta[[ yCatBra ]] ) ) ) xBetaBra <- colSums( allCoefBra * mXValBra ) pfunBra <- exp( xBetaBra[ yCatBra ] + lambda[ yCatBra ] * IV[ yCatBra ] )/ ( sum( exp( xBetaBra + lambda * IV ) ) ) semEla <- ( xCoef[1] + 2 * xCoef[2] * xVal[ yCat ] ) * xVal[ yCat ] * ( pfunBra * ( pfun - pfun^2 ) * 1/lambda[ yCatBra ] + pfun^2 * ( pfunBra - pfunBra^2 ) * lambda[ yCatBra ] * IV[ yCatBra ] ) } if( method == "binary" || method == "MNL" ){ derivCoef <- rep( 0, length( allCoef ) ) derivCoef[ xPos[1] ] <- dfun * xVal if( length( xPos ) == 2 ) { derivCoef[ xPos[2] ] <- dfun * 2 * xVal^2 } } else if( method == "CondL" ){ derivCoef <- rep( 0, length( allCoef ) ) derivCoef[ xPos[1] ] <- ( pfun - pfun^2 ) * xVal[ yCat ] if( length( xPos ) == 2 ) { derivCoef[ xPos[2] ] <- ( pfun - pfun^2 ) * 2 * xVal[ yCat ]^2 } } else{ derivCoef <- rep( 0, length( allCoef ) ) derivCoef[ xPos[1] ] <- ( pfunBra * ( pfun - pfun^2 ) / lambda[ yCatBra ] + pfun^2 * ( pfunBra - pfunBra^2 ) * lambda[ yCatBra ] * IV[ yCatBra ] ) * xVal[ yCat ] if( length( xPos ) == 2 ) { derivCoef[ xPos[2] ] <- ( pfunBra * ( pfun - pfun^2 ) / lambda[ yCatBra ] + pfun^2 * ( pfunBra - pfunBra^2 ) * lambda[ yCatBra ] * IV[ yCatBra ] ) * 2 * xVal[ yCat ]^2 } } vcovCoef <- diag( allCoefSE^2 ) semElaSE <- drop( sqrt( t( derivCoef ) %*% vcovCoef %*% derivCoef ) ) result <- c( semEla = semEla, stdEr = semElaSE ) return( result ) } @ with \begin{itemize} \item arguments \code{allCoef} $= (\beta_0, \ldots, \beta_K )^\top$ a vector of all coefficients from a binary logit model, \code{allCoef} $= ( \gamma_0, \ldots, \gamma_T )^\top$ a vector of all coefficients from a conditional logit model, \code{allCoefBra} $= ( \beta_0, \ldots, \beta_K )^\top$, a vector of all coefficients from a nested logit model, where \code{allCoefBra}, the coefficients at the branch level, must always be included in combination with \code{allCoef}, the coefficients at the twig level, or \code{allCoef} $= ( \beta_{01}, \ldots, \beta_{K1}, \ldots, \beta_{0P}, \ldots, \beta_{KP} )$ a vector of the $P$ sets of coefficients from the multi-nomial logit regression which does \textbf{not} include any values for the reference category; \item argument \code{allXVal} $= ( 1, \ldots, \bar{x}_K )^\top$, a vector of all corresponding sample means and 1 for the intercept for the binary or multi-nomial logit, \code{allXVal} $= ( 1, \ldots, \bar{z}_{T,1}, \ldots, 1, \ldots, \bar{z}_{T,P})$ a vector of the $P$ sets of explanatory variables from the conditional logit, or \code{allXVal} $=((( 1, \ldots, \bar{x}_{K,p,o} ), \ldots, (1, \ldots, \bar{x}_{K,P,o}))^\top, \ldots, (( 1, \ldots, \bar{x}_{K,p,O}), \ldots, ( 1, \ldots, \bar{x}_{K,P,O}))^\top )^\top$, a list where the list elements are the corresponding matrices of the sample means for each nest at the twig level,% \footnote{The reason why we make an exception from the usual vector in the case of the nested logit model is that nests, $o$, can have different dimensions wrt $P$.} and which must always be combined with the corresponding values at the branch level in \code{allXVal}; \item argument \code{allCoefSE} $= (\se(\beta_0), \ldots, \se( \beta_K))^\top$, a vector of standard errors for all coefficients of a binary or conditional logit regression or of the standard errors of the $P$ sets of coefficients in a multi-nomial logit regression; \item argument \code{xPos} $= k+1$ or $(k+1, k+2)^\top$ a scalar or vector of two elements indicating the position of the coefficients of interest in vectors \code{allCoef} and \code{allXVal}; \item argument \code{method = c("binary", "MNL", "CondL", "NestL" )} indicating the estimator used, where option \code{"CondL"} can also be used to calculate the semi-elasticity and standard error of a nested logit at the branch level and where option \code{"NestL"} estimates the semi-elasticity and standard error of the combined probabilities at the branch and twig level; \item argument \code{"lambda"} is a vector of the $O$ parameter of the correlation coefficient between the the branch specific alternatives of interest, it is only relevant under option \code{"NestL"} and should be set to $1$ in case the estimation coefficients in \code{allCoef} are already corrected by $\frac{ \beta_k }{ \lambda_o }$; \item finally, \code{yCat} a scalar identifying which of the $P$ output categories is of interest, only in combination with \code{method = "MNL"}, \code{"CondL"}, and \code{"NestL"} (for the twig level), and \code{yCatBra} at the branch level. \end{itemize} << >>= # Example ela6a <- logEla( allCoef = c( 0.445, 0.03, 0.00002, 0.067, 0.89, 0.124 ), allXVal = c( 1, 3.3, 4.5, 2.34, 0.1, 0.987 ), xPos = 2 ) ela6a ela6b <- logEla( allCoef = c( 0.445, 0.03, 0.00002, 0.067, 0.89, 0.124 ), allXVal = c( 1, 3.3, 4.5, 2.24, 0.1, 0.987 ), allCoefSE = c( 0.001, 0.02, 0.000002, 0.05, 1.2, 0.03 ), xPos = 2 ) ela6b # Example ela7a <- logEla( allCoef = c( 0.445, 0.03, 0.00002, 0.067, 0.89, 0.124 ), allXVal = c( 1, 3.3, 3.3^2, 2.34, 0.1, 0.987 ), xPos = c( 2, 3 ) ) ela7a ela7b <- logEla( allCoef = c( 0.445, 0.03, 0.00002, 0.067, 0.89, 0.124 ), allXVal = c( 1, 3.3, 3.3^2, 2.34, 0.1, 0.987 ), allCoefSE = c( 0.001, 0.02, 0.000002, 0.05, 1.2, 0.03 ), xPos = c( 2, 3 ) ) ela7b # Example ela8a <- logEla( allCoef = c( 0.2, 0.3, 0.5, -0.2, 0.03, 0.6 ), allXVal = c( 1, 8.4, 0.06 ), xPos = 3, method = "MNL", yCat = 2 ) ela8a ela8b <- logEla( allCoef = c( 0.2, 0.3, 0.5, -0.2, 0.03, 0.6 ), allXVal = c( 1, 8.4, 0.06 ), allCoefSE = c( 0.002, 0.003, 0.004, 0.006, 0.00001, 0.08 ), xPos = 3, method = "MNL", yCat = 2 ) ela8b # Example ela9a <- logEla( allCoef = c( 0.2, 0.3, 0.5, -0.2, 0.03, 0.6 ), allXVal = c( 1, 0.04, 0.0016 ), xPos = c( 2, 3 ), method = "MNL", yCat = 2 ) ela9a ela9b <- logEla( allCoef = c( 0.2, 0.3, 0.5, -0.2, 0.03, 0.6 ), allXVal = c( 1, 0.04, 0.0016 ), allCoefSE = c( 0.002, 0.003, 0.004, 0.006, 0.00001, 0.08 ), xPos = c( 2, 3 ), method = "MNL", yCat = 2 ) ela9b # Example ela10a <- logEla( allCoef = c( 0.445, 0.03, 0.00002 ), allXVal = c( 1, 3.3, 4.5, 1, 0.1, 0.987 ), xPos = 2, method = "CondL", yCat = 2 ) ela10a ela10b <- logEla( allCoef = c( 0.445, 0.03, -0.002 ), allXVal = c( 1, 0.3, 0.09, 1, 0.1, 0.01 ), xPos = c( 2, 3 ), method = "CondL", yCat = 2 ) ela10b # Example ela11a <- logEla( allCoef = c( 0.445, 0.03, 0.00002 ), allXVal = c( 1, 3.3, 4.5, 1, 0.1, 0.987 ), allCoefSE = c( 0.002, 0.003, 0.004 ), xPos = 2, method = "CondL", yCat = 2 ) ela11a ela11b <- logEla( allCoef = c( 0.445, 0.03, -0.002 ), allXVal = c( 1, 0.3, 0.09, 1, 0.1, 0.01 ), allCoefSE = c( 0.002, 0.003, 0.004 ), xPos = c( 2, 3 ), method = "CondL", yCat = 2 ) ela11b # Example matrix1 <- matrix( c( 1, 2.5, 0.3, 0.09, 1, 0.33, 0.9, 1.8 ), nrow = 4 ) matrix2 <- matrix( c( 1, 2.8, 0.099, 0.211 ), nrow = 4 ) ela12a <- logEla( allCoefBra = c( 0.445, 0.03, -0.002 ), allCoef = c( 1.8, 0.005, -0.12, 0.8 ), allXValBra = c( 1, 3.3, 4.5, 1, 0.1, 0.987 ), allXVal = list( matrix1, matrix2 ), xPos = 2, yCatBra = 1, yCat = 2, lambda = c( 0.8, 1 ), method = "NestedL" ) ela12a matrix1 <- matrix( c( 1, 0.3, 0.09, 0.09, 1, 0.33, 0.1089, 1.8 ), nrow = 4 ) matrix2 <- matrix( c( 1, 0.31, 0.099, 0.211 ), nrow = 4 ) ela12b <- logEla( allCoefBra = c( 0.445, 0.03, -0.002 ), allCoef = c( 1.8, 0.005, -0.12, 0.8 ), allXValBra = c( 1, 3.3, 4.5, 1, 0.1, 0.987 ), allXVal = list( matrix1, matrix2 ), xPos = 2, yCatBra = 1, yCat = 2, lambda = c( 0.8, 1 ), method = "NestedL" ) ela12b # Example matrix1 <- matrix( c( 1, 2.5, 0.3, 0.09, 1, 0.33, 0.9, 1.8 ), nrow = 4 ) matrix2 <- matrix( c( 1, 2.8, 0.099, 0.211 ), nrow = 4 ) ela13a <- logEla( allCoefBra = c( 0.445, 0.03, -0.002 ), allCoef = c( 1.8, 0.005, -0.12, 0.8 ), allXValBra = c( 1, 3.3, 4.5, 1, 0.1, 0.987 ), allXVal = list( matrix1, matrix2 ), allCoefSE = c( 0.001, 0.089, 0.0003, 0.12 ), xPos = 2, yCatBra = 1, yCat = 2, lambda = c( 0.8, 1 ), method = "NestedL" ) ela13a matrix1 <- matrix( c( 1, 0.3, 0.09, 0.09, 1, 0.33, 0.1089, 1.8 ), nrow = 4 ) matrix2 <- matrix( c( 1, 0.31, 0.099, 0.211 ), nrow = 4 ) ela13b <- logEla( allCoefBra = c( 0.445, 0.03, -0.002 ), allCoef = c( 1.8, 0.005, -0.12, 0.8 ), allXValBra = c( 1, 3.3, 4.5, 1, 0.1, 0.987 ), allXVal = list( matrix1, matrix2 ), allCoefSE = c( 0.001, 0.089, 0.0003, 0.12 ), xPos = 2, yCatBra = 1, yCat = 2, lambda = c( 0.8, 1 ), method = "NestedL" ) ela13b @ %============================================================================================= \subsection{Semi-elasticities from interval-coded explanatory variables} Logit model, where the $k$th explanatory variable is interval-coded: \begin{equation} \Pr( y = 1 | x ) = \frac{ \exp{( \beta_0 + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j x_j + \sum_{m \in \{1, \ldots, M \} \setminus m^* } \delta_m D_m )}} { 1 + \exp{( \beta_0 + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j x_j + \sum_{m \in \{1, \ldots, M \} \setminus m^* } \delta_m D_m )}}, \label{eq:inLogit} \end{equation} with \begin{equation} D_m = \begin{cases} 1 & \text{ if } b_{m-1} < x_k \leq b_m \\ 0 & \text{ otherwise} \end{cases} \quad \forall \; m = 1, \ldots , M, \end{equation} The semi-elasticities from an interval coded variable in the case of a binary logit function can be calculated as in section~\ref{sec:probInt}, where equation~(\ref{eq:probElaInt}) modifies to \begin{equation} \epsilon_{km} \approx 2 \left( \frac{ \exp_{m+1}{(\cdot)}}{1 + \exp_{m+1}{(\cdot)}} - \frac{ \exp_m{(\cdot)}}{1 + \exp_m{(\cdot)}} \right) \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \label{eq:logSem} \end{equation} where \begin{equation} \exp_{m+1}{(\cdot)} \equiv \exp \left( \beta_0 + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j x_j + \delta_{m+1} \right) \end{equation} and \begin{equation} \exp_m{(\cdot)} \equiv \exp \left( \beta_0 + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j x_j + \delta_m \right). \end{equation} For the multi-nomial logit, equation~(\ref{eq:logSem}) modifies to \begin{equation} \epsilon_{km,p} \approx 2 \left( \frac{ \exp_{m+1, p}{(\cdot)}} {1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp_{m+1, p}{(\cdot)}} - \frac{ \exp_{m, p}{(\cdot)}} {1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp_{m,p}{(\cdot)}} \right) \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \label{eq:MNlogSem} \end{equation} where $\exp_{m, p}$ is defined as \begin{equation} \exp_{m,p}{(\cdot)} \equiv \exp \left( \beta_{0,p} + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_{j,p} x_j + \delta_{m,p} \right). \end{equation} and $\exp_{m+1, p}$ accordingly. For the conditional logit, equation~(\ref{eq:logSem}) modifies to \begin{equation} \epsilon_{tm,p} \approx 2 \left( \frac{ \exp_{m+1, p}{(\cdot)}}{\sum_{p \in C } \exp_{m+1, p}{(\cdot)}} - \frac{ \exp_{m, p}{(\cdot)}}{\sum_{p \in C } \exp_{m,p}{(\cdot)}} \right) \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \label{eq:CondlogSem} \end{equation} where $\exp_{m, p}$ is defined as \begin{equation} \exp_{m,p}{(\cdot)} \equiv \exp \left( \gamma_0 + \sum_{j \in \{1, \ldots, T \} \setminus t } \gamma_j z_{j,p} + \delta_m \right). \end{equation} and $\exp_{m+1, p}$ accordingly. An approximate standard error of the semi-elasticity of interval-coded explanatory variables of binary logit models can be obtained by using the Delta-method: \begin{align} \se( \epsilon_{km}) & = \sqrt{ \frac{ \partial \epsilon_{km} }{ \partial \left( \beta^{\top} \; \delta^{\top} \right) } \; \var \left( \begin{array}{c} \beta \\ \delta \end{array} \right) \; \frac{ \partial \epsilon_{km} }{ \partial \left( \begin{array}{c} \beta \\ \delta \end{array} \right) } }, \label{eq:logElaIntSE} \end{align} where $\se( \epsilon_{km} )$ indicates the (approximate) standard error of $\epsilon_{km}$, $\var \left( \left( \beta^{\top} \; \delta^{\top} \right)^{\top} \right)$ indicates the variance-covariance matrix of the estimates of the coefficient vector~% $\left( \beta^{\top} \; \delta^{\top} \right)^{\top}$, and the gradient vector is: \begin{equation} \frac{ \partial \epsilon_{km}} { \partial \begin{pmatrix} \beta \\ \delta \end{pmatrix}} = \sum_{m = 1}^{M-1} w_m \frac{ \partial \epsilon_{km}} { \partial \begin{pmatrix} \beta \\ \delta \end{pmatrix}} \end{equation} with the elements of the gradient vector~% $\partial \epsilon_{km} / \partial ( \beta^{\top} \delta^{\top} )^{\top}$: \begin{align} \frac{ \partial \epsilon_{km} }{ \partial \beta_0 } & = 2 \left( \frac{ \exp_{m+1}(\cdot)}{ [ 1 + \exp_{m+1}{( \cdot )} ]^2 } - \frac{ \exp_m (\cdot)}{ [ 1 + \exp_m{( \cdot )} ]^2 } \right) \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \\[1em] \frac{ \partial \epsilon_{km} }{ \partial \beta_j } & = 2 \left( \frac{ \exp_{m+1}(\cdot) x_j}{ [ 1 + \exp_{m+1}{( \cdot )} ]^2 } - \frac{ \exp_m (\cdot) x_j}{ [ 1 + \exp_m{( \cdot )} ]^2 } \right) \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \; \forall \; j \in \{ 1, \ldots, K \} \setminus k \\[1em] \frac{ \partial \epsilon_{km} }{ \partial \delta_m } & = -2 \frac{ \exp_m (\cdot)}{[ 1 + \exp_m (\cdot )]^2 } \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \\[1em] \frac{ \partial \epsilon_{km} }{ \partial \delta_{m+1} } & = 2 \frac{ \exp_{m+1} (\cdot)}{[ 1 + \exp_{m+1} (\cdot )]^2 } \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \\[1em] \frac{ \partial \epsilon_{km} }{ \partial \delta_n } & = 0 \; \forall \; n \in \{ 1, \ldots , M \} \setminus \{ m, m+1 \} \label{eq:logGrad} \end{align} In the case of the multi-nomial logit the elements of the gradient vector~% $\partial \epsilon_{km,p} / \partial ( \beta^{\top} \delta^{\top} )^{\top}$ modify to \begin{align} \frac{ \partial \epsilon_{km,p}}{\partial \beta_{0,p}} =& \, 2 \left( \frac{ \exp_{m + 1, p}( \cdot ) \cdot ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ m + 1, o } ( \cdot ) ) } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m+1, p }( \cdot ) ]^2 } \right. \\ & \left. - \frac{ \exp_{m, p}( \cdot ) \cdot ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ m, o } ( \cdot )) } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \right) \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \nonumber \label{eq:MNLgrad1} \\[1em] \frac{ \partial \epsilon_{km,p}}{\partial \beta_{0,o}} =& \, 2 \left( \frac{ \exp_{m,p}( \cdot ) \exp_{m,o}( \cdot )} {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } - \frac{ \exp_{m+1,p}( \cdot ) \exp_{m+1,o}( \cdot )} {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m+1, p }( \cdot ) ]^2 } \right) \\ & \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \;\;\; \forall \; o \neq p \nonumber \\[1em] \frac{ \partial \epsilon_{km,p}}{\partial \beta_{j,p}} =& \, 2 \left( \frac{ \exp_{m + 1, p}( \cdot ) \cdot x_j ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ m + 1, o } ( \cdot ) ) } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m+1, p }( \cdot ) ]^2 } \right. \\ & \left. - \frac{ \exp_{m, p}( \cdot ) \cdot x_j ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ m, o } ( \cdot )) } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \right) \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \;\;\; \forall \; j \in \{ 1, \ldots, K \} \setminus k \nonumber \\[1em] \frac{ \partial \epsilon_{km,p}}{\partial \beta_{j,o}} =& \, 2 \left( \frac{ \exp_{m,p}( \cdot ) \exp_{m,o}( \cdot ) \cdot x_j } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } - \frac{ \exp_{m+1,p}( \cdot ) \exp_{m+1,o}( \cdot ) \cdot x_j } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m+1, p }( \cdot ) ]^2 } \right) \\ & \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \;\;\; \forall \; o \neq p, \; j \in \{ 1, \ldots, K \} \setminus k \nonumber \\[1em] \frac{ \partial \epsilon_{km}}{\partial \delta_{m,p}} =& - 2 \frac{ \exp_{m, p}( \cdot ) \cdot ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ m, o } ( \cdot )) } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \\[1em] \frac{ \partial \epsilon_{km,p}}{\partial \delta_{m,o}} =& \, 2 \frac{ \exp_{m,p}( \cdot ) \exp_{m,o}( \cdot )} {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \;\;\; \forall \; o \neq p \\[1em] \frac{ \partial \epsilon_{km,p}}{\partial \delta_{m+1,p}} =& \, 2 \frac{ \exp_{m+1, p}( \cdot ) \cdot ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ m + 1, o } ( \cdot )) } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m + 1, p }( \cdot ) ]^2 } \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \\[1em] \frac{ \partial \epsilon_{km,p}}{\partial \delta_{m + 1,o}} =& -2 \frac{ \exp_{m + 1,p}( \cdot ) \exp_{m + 1,o}( \cdot )} {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m + 1, p }( \cdot ) ]^2 } \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \;\;\; \forall \; o \neq p \\[1em] \frac{ \partial \epsilon_{km,p} }{ \partial \delta_{n,p}} =& \; 0 \;\; \forall \; n \in \{ 1, \ldots , M \} \setminus \{ m, m+1 \} \\[1em] \frac{ \partial \epsilon_{km,p} }{ \partial \delta_{n,o}} =& \; 0 \;\; \forall \; n \in \{ 1, \ldots , M \} \setminus \{ m, m+1 \}, \; o \neq p \label{eq:MNLgradn} \end{align} In the case of the conditional logit the elements of the gradient vector~% $\partial \epsilon_{tm,p} / \partial ( \gamma^{\top} \delta^{\top} )^{\top}$ modify to \begin{align} \frac{ \partial \epsilon_{tm} }{ \partial \gamma_0 } =& \; 2 \left( \frac{ \exp_{m+1,p}(\cdot) \cdot \sum_{p\in C} \exp_{m+1,p}(\cdot) - \exp_{m+1,p}(\cdot) \cdot \sum_{p\in C} \exp_{m+1,p}(\cdot) } { \left( \sum_{p\in C} \exp_{m+1,p}(\cdot) \right)^2 } - \right. \label{eq:condLogGrad1}\\ & \left. - \frac{ \exp_{m,p}(\cdot) \cdot \sum_{p\in C} \exp_{m,p}(\cdot) - \exp_{m,p}(\cdot) \cdot \sum_{p\in C} \exp_{m,p}(\cdot) } { \left( \sum_{p\in C} \exp_{m,p}(\cdot) \right)^2 } \right) \cdot \frac{b_m}{b_{m+1} - b_{m-1}} = 0 \nonumber \\[1em] \frac{ \partial \epsilon_{tm} }{ \partial \gamma_j } =& \; 2 \left( \frac{ \exp_{m+1,p}(\cdot) z_{j,p} \cdot \sum_{p\in C} \exp_{m+1,p}(\cdot) - \exp_{m+1,p}(\cdot) \cdot \sum_{p\in C} \exp_{m+1,p}(\cdot) z_{j,p} } { \left( \sum_{p\in C} \exp_{m+1,p}(\cdot) \right)^2 } \right. \\ & \left. - \frac{ \exp_{m,p}(\cdot) z_{j,p} \cdot \sum_{p\in C} \exp_{m,p}(\cdot) - \exp_{m,p}(\cdot) \cdot \sum_{p\in C} \exp_{m,p}(\cdot) z_{j,p} } { \left( \sum_{p\in C} \exp_{m,p}(\cdot) \right)^2 } \right) \nonumber \\[1em] & \cdot \frac{b_m}{b_{m+1} - b_{m-1}} \; \forall \; j \in \{ 1, \ldots, T \} \setminus t \nonumber \\[1em] \frac{ \partial \epsilon_{tm} }{ \partial \delta_m } =& \; -2 \frac{ \exp_{m,p}(\cdot) \cdot \sum_{p\in C} \exp_{m,p}(\cdot) - \exp_{m,p}(\cdot) \cdot \sum_{p\in C} \exp_{m,p}(\cdot) } { \left( \sum_{p\in C} \exp_{m,p}(\cdot) \right)^2 } \\ & \cdot \frac{b_m}{b_{m+1} - b_{m-1}} = 0 \nonumber \\[1em] \frac{ \partial \epsilon_{tm} }{ \partial \delta_{m+1} } =& \; 2 \frac{ \exp_{m+1,p}(\cdot) \cdot \sum_{p\in C} \exp_{m+1,p}(\cdot) - \exp_{m+1,p}(\cdot) \cdot \sum_{p\in C} \exp_{m+1,p}(\cdot) } { \left( \sum_{p\in C} \exp_{m+1,p}(\cdot) \right)^2 } \\ & \cdot \frac{b_m}{b_{m+1} - b_{m-1}} = 0 \nonumber \\[1em] \frac{ \partial \epsilon_{km} }{ \partial \delta_n } =& \; 0 \; \forall \; n \in \{ 1, \ldots , M \} \setminus \{ m, m+1 \} \label{eq:condLogGrad} \end{align} Using the helper functions \code{elaIntWeights} (equation~\ref{eq:linElaIntWeights}), \code{elaIntBounds} (equation~\ref{eq:linElaIntBM}), \code{checkXPos} (checking argument \code{xPos}, see appendix~\ref{sec:helperFunctions}), and \code{checkXBeta} (checking $\beta_0 + \sum_{j\in\{1, \ldots, K \} \setminus k} \beta_j x_j + \sum_{m \in \{1,...,M \} \setminus m^*} \delta_m D_m$ for plausible values, see appendix~\ref{sec:helperFunctions}), and function \code{linElaInt} (section~\ref{sec:linInterval}), the following code defines a function that calculates the semi-elasticity defined in equation~(\ref{eq:logSem}) and its approximate standard error defined in equations~(\ref{eq:logElaIntSE}) to~(\ref{eq:logGrad}) for the binary logit function, equations~(\ref{eq:MNlogSem}) and (\ref{eq:MNLgrad1}) to (\ref{eq:MNLgradn}) for the multi-nomial logit function, and equations~(\ref{eq:CondlogSem}) and (\ref{eq:condLogGrad1}) to (\ref{eq:condLogGrad}) for the conditional logit, respectively: << >>= logElaInt <- function( allCoef, allXVal, xPos, xBound, yCat = NA, allCoefSE = rep( NA, length( allCoef ) ), method = "binary" ){ # number of coefficients if( method == "binary" || method == "MNL" ){ mCoef <- matrix( allCoef, nrow = length( allXVal )) nCoef <- dim( mCoef )[1] pCoef <- dim( mCoef )[2] # checking arguments if( length( allXVal ) != nCoef ) { stop( "arguments 'allCoef' and 'allXVal' must have the same length" ) } } else{ nCoef <- length( allCoef ) mXVal <- matrix( allXVal, nrow = nCoef ) pCoef <- dim( mXVal )[2] # checking arguments if( dim( mXVal )[1] != nCoef ) { stop( "arguments 'allCoef' and 'allXVal' must have the same length" ) } } # number of intervals nInt <- length( xPos ) checkXPos( xPos, minLength = 2, maxLength = nCoef, minVal = 0, maxVal = nCoef, requiredVal = 0 ) if( method == "binary" || method == "MNL" ){ if( any( allXVal[ xPos ] < 0 ) ) { stop( "all elements of argument 'allXVal'", " that are indicated by argument 'xPos'", " (i.e., the shares of observations in each interval)", " must be non-negative" ) } if( sum( allXVal[ xPos ] > 1 ) ) { stop( "the sum of the elements of argument 'allXVal'", " that are indicated by argument 'xPos'", " (i.e., the shares of observations in each interval)", " must not be larger than one" ) } } else{ for( p in 1:pCoef ){ if( any( mXVal[ xPos, p ] < 0 ) ) { stop( "all elements of argument 'allXVal'", " that are indicated by argument 'xPos'", " (i.e., the shares of observations in each interval)", " must be non-negative" ) } if( sum( mXVal[ xPos, p ] > 1 ) ) { stop( "the sum of the elements of argument 'allXVal'", " that are indicated by argument 'xPos'", " (i.e., the shares of observations in each interval)", " must not be larger than one" ) } } } # check 'xBound' and replace infinite values xBound <- elaIntBounds( xBound, nInt ) # vector of probabilities of y=1 for each interval and # vector of shares of observations in each interval xBeta <- matrix( rep( rep( NA, nInt ), pCoef ), ncol = pCoef ) if( method == "binary" || method == "MNL" ){ shareVec <- rep( NA, nInt ) for( p in 1:pCoef ){ for( i in 1:nInt ){ allXValTemp <- replace( allXVal, xPos, 0 ) if( xPos[i] != 0 ) { allXValTemp[ xPos[i] ] <- 1 shareVec[i] <- allXVal[ xPos[i] ] } xBeta[i,p] <- sum( mCoef[ ,p] * allXValTemp ) } } shareVec[ xPos == 0 ] <- 1 - sum( shareVec[ xPos != 0 ] ) } else{ shareVec <- matrix( rep( rep( NA, nInt ), pCoef ), ncol = pCoef ) for( p in 1:pCoef ){ for( i in 1:nInt ){ allXValTemp <- replace( mXVal[ ,p], xPos, 0 ) if( xPos[i] != 0 ) { allXValTemp[ xPos[i] ] <- 1 shareVec[i,p] <- mXVal[ xPos[i], p ] } xBeta[i,p] <- sum( allCoef * allXValTemp ) } shareVec[ xPos == 0, p ] <- 1 - sum( shareVec[ xPos != 0, p ] ) } shareVec <- shareVec[ , yCat ] } #checkXBeta( xBeta ) #Please check this one with a matrix if( method == "binary" ){ expVec <- as.vector( exp( xBeta )/( 1 + exp( xBeta ) ) ) } else if( method == "MNL" ){ expVec <- as.vector( exp( xBeta[ , yCat ])/( 1 + rowSums( exp( xBeta ) ) ) ) } else{ expVec <- as.vector( exp( xBeta[ , yCat ])/( rowSums( exp( xBeta ) ) ) ) } # weights weights <- elaIntWeights( shareVec ) # calculation of the semi-elasticity semEla <- linElaInt( expVec, shareVec, xBound ) ### calculation of its standard error # partial derivatives of each semi-elasticity around each boundary # w.r.t. all estimated coefficients if( method == "binary" ){ gradM <- matrix( 0, nCoef, nInt - 1 ) gradExpVec <- exp( xBeta )/( 1 + exp( xBeta ) )^2 for( m in 1:( nInt - 1 ) ) { gradM[ -xPos, m ] <- 2 * ( gradExpVec[m+1] - gradExpVec[m] ) * allXVal[ -xPos ] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) gradM[ xPos[m], m ] <- - 2 * gradExpVec[m] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) gradM[ xPos[m+1], m ] <- 2 * gradExpVec[m+1] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) } } else if( method == "MNL" ){ gradM <- array( 0, c( nCoef, nInt - 1, pCoef ) ) gradExpVecP <- ( exp( xBeta[ , yCat ] ) * ( 1 + rowSums( exp( xBeta[ , -yCat, drop = FALSE ] ) ) ) )/ ( 1 + rowSums( exp( xBeta ) ) )^2 for( p in 1:pCoef ){ gradExpVecO <- ( exp( xBeta[ , yCat ] ) * exp( xBeta[ , p] ) )/ ( 1 + rowSums( exp( xBeta ) ) )^2 for( m in 1:( nInt - 1 ) ) { if( p == yCat ){ gradM[ -xPos, m, p ] <- 2 * ( gradExpVecP[m+1] - gradExpVecP[m] ) * allXVal[ -xPos ] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) gradM[ xPos[ m ], m, p ] <- - 2 * gradExpVecP[m] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) gradM[ xPos[ m + 1 ], m, p ] <- 2 * gradExpVecP[m+1] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) } else { gradM[ -xPos, m, p ] <- 2 * ( gradExpVecO[m] - gradExpVecO[m+1] ) * allXVal[ -xPos ] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) gradM[ xPos[ m ], m, p ] <- 2 * gradExpVecO[m] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) gradM[ xPos[ m + 1 ], m, p ] <- - 2 * gradExpVecO[m+1] * xBound[m+1] / ( xBound[m+2] - xBound[m] ) } } } gradM <- apply( gradM, 2, function( x ) x ) } else{ gradM <- matrix( 0, nCoef, nInt - 1 ) for( m in 1:( nInt - 1 ) ) { gradM[ -xPos, m ] <- 2 * ( ( exp( xBeta[ m+1, yCat ] ) * mXVal[ -xPos, yCat ] * sum( exp( xBeta[ m+1, ] ) ) - exp( xBeta[ m+1, yCat ] ) * rowSums( exp( xBeta[ m+1, ] ) * mXVal[ -xPos, , drop = FALSE ] ) )/ ( sum( exp( xBeta[ m+1, ] ) ) )^2 - ( exp( xBeta[ m, yCat ] ) * mXVal[ -xPos, yCat ] * sum( exp( xBeta[ m, ] ) ) - exp( xBeta[ m, yCat ] ) * rowSums( exp( xBeta[ m, ] ) * mXVal[ -xPos, , drop = FALSE ] ) )/ ( sum( exp( xBeta[ m, ] ) ) )^2 ) * xBound[m+1] / ( xBound[m+2] - xBound[m] ) gradM[ xPos[m], m ] <- 0 gradM[ xPos[m+1], m ] <- 0 } } # partial derivative of the semi-elasticity # w.r.t. all estimated coefficients derivCoef <- rep( 0, length( allCoef ) ) for( m in 1:( nInt - 1 ) ){ derivCoef <- derivCoef + weights[m] * gradM[,m] } # variance-covariance matrix of the coefficiencts vcovCoef <- diag( allCoefSE^2 ) # standard error of the (average) semi-elasticity semElaSE <- drop( sqrt( t( derivCoef ) %*% vcovCoef %*% derivCoef ) ) # prepare object that will be returned result <- c( semEla[1], stdEr = semElaSE ) return( result ) } @ where argument \code{allCoef} $= ( \beta_0, \ldots, \beta_{k-1}, \delta_1, \ldots, \delta_{m^*-1}, \delta_{m^*+1}, \ldots, \delta_M, \beta_{k+1}, \ldots, \beta_K )^{\top}$, is a vector of all coefficients from the binary or conditional logit regression or \code{allCoef} $= ( \beta_{0,1}, \ldots, \beta_{k-1,1}, \delta_{1,1}, \ldots, \delta_{m^*-1,1}, \delta_{m^*+1,1}, \ldots, \delta_{M,1}, \beta_{k+1,1}, \ldots, \beta_{K,1}, \ldots, \\ \beta_{0,P}, \ldots, \beta_{k-1,P}, \delta_{1,P}, \ldots, \delta_{m^*-1,P}, \delta_{m^*+1,P}, \ldots, \delta_{M,P}, \beta_{k+1,P}, \ldots, \beta_{K,P} )^\top$ a vector of all the $P$ sets of coefficients from the multi-nomial logit regression which does \textbf{not} include any values for the reference category; argument \code{allXVal} $= ( 1, x_1, \ldots, x_{k-1},$ $s_1, \ldots, s_{m^*-1}, s_{m^*+1}, \ldots, s_M,$ $x_{k+1}, \ldots, x_K )^{\top}$ is a vector of corresponding values of the explanatory variables (including shares of observations in each interval of variable~$x_k$ except for the reference interval~$m^*$) or \code{allXVal} $= ( 1, z_{1,1}, \ldots, z_{1,T}, \ldots, 1, z_{1,P}, \ldots, z_{T,P})$ is a vector of the $P$ sets of explanatory variables from the conditional logit (including shares of observations in each interval of variable~$z_{t,p}$ except for the reference interval~$m^*$); argument \code{xPos} $= ( k+1, \ldots, k + m^* - 1, 0, k + m^* , \ldots , k + M - 1 )$ is a vector indicating the positions of the coefficients~$\delta_1, \ldots, \delta_M$ and shares~$s_1, \ldots, s_M$ of each interval in the vectors \code{allCoef} and \code{allXVal}, whereas the position of the reference interval~$m^*$ is set to zero; argument \code{allCoefSE} is a vector of the standard errors of all coefficients in \code{allCoef} or of the $P$ sets of the standard errors of the coefficients from the multi-nomial logit regression which does \textbf{not} include any values for the reference category; argument \code{method = "binary", "MNL", "CondL"} identifies the estimator chosen, where \code{"binary"} indicates the binary logit estimator, \code{"MNL"} indicates the multi-nomial logit estimator, and \code{"CondL"} indicates the conditional logit estimator; argument \code{yCat}, only in combination with option \code{"MNL"} or \code{"CondL"}, indicates the $p$th output category of interest; and argument \code{xBound} $= ( b_0 , \ldots , b_M )^{\top}$ indicates the boundaries of the intervals as in function \code{linElaInt}. << >>= # Example ela8a <- logElaInt( allCoef = c( 0.33, 0.22, 0.05, 0.6 ), allXVal = c( 1, 0.4, 0.12, 0.13 ), xPos = c( 2, 0, 3, 4 ), xBound = c( 0, 500, 1000, 1500, Inf ) ) ela8a ela8b <- logElaInt( allCoef = c( 0.33, 0.22, 0.05, 0.6 ), allXVal = c( 1, 0.4, 0.12, 0.13 ), xPos = c( 2, 0, 3, 4 ), xBound = c( 0, 500, 1000, 1500, Inf ), allCoefSE = c( 0.003, 0.045, 0.007, 0.009 ) ) ela8b # Example ela9a <- logElaInt( allCoef = c( 0.2, 0.3, 0.5, -0.2, 0.03, 0.6 ), allXVal = c( 1, 0.4, 0.12 ), xPos = c( 2, 0, 3 ), xBound = c( 0, 500, 1000, Inf ), yCat = 2, method = "MNL" ) ela9a ela9b <- logElaInt( allCoef = c( 0.2, 0.3, 0.5, -0.2, 0.03, 0.6 ), allXVal = c( 1, 0.4, 0.12 ), xPos = c( 2, 0, 3 ), xBound = c( 0, 500, 1000, Inf ), yCat = 2, allCoefSE = c( 0.003, 0.045, 0.007, 0.009, 0.0008, 0.9 ), method = "MNL" ) ela9b # Example ela10a <- logElaInt( allCoef = c( 1.33, 0.022, 0.58, 1.6 ), allXVal = c( 1, 0.4, 0.12, 0.0002, 1, 0.28, 0.01, 0.000013 ), xPos = c( 2, 0, 3 ), xBound = c( 0, 1000, 1500, Inf ), yCat = 2, method = "CondL" ) ela10a ela10b <- logElaInt( allCoef = c( 1.33, 0.022, 0.58, 1.6 ), allXVal = c( 1, 0.4, 0.12, 0.0002, 1, 0.28, 0.01, 0.000013 ), xPos = c( 2, 0, 3 ), xBound = c( 0, 1000, 1500, Inf ), allCoefSE = c( 0.003, 0.045, 0.007, 0.009 ), yCat = 2, method = "CondL" ) ela10b @ %============================================================================================= \subsection{Effects of continuous variables when they change between discrete intervals} As in section~\ref{seq:logit}, we assume the following model \begin{equation} \Pr( y = 1 | x ) = \frac{ \exp{( \beta_0 + \sum_{k=1}^K \beta_k x_k )}} { 1 + \exp{( \beta_0 + \sum_{k=1}^K \beta_k x_k )}}, \end{equation} where $y \in \{0, 1 \}$ is a binary dependent variable, $x = (x_1, \ldots, x_K )^{\top}$ is a vector of $K$ continuous explanatory variables, and $\beta = ( \beta_0, \ldots, \beta_K )^{\top}$ is a vector of $K + 1$ unknown coefficients. In order to derive the (approximate) effect of a continuous variable $x_k$ on $Pr(y = 1)$ given that $x_k$ changes between $M \geq 2$ intervals, we modify equation~(\ref{eq:linEffect}) into \begin{align} E_{k,ml} = & E[ y | b_{m-1} < x_k \leq b_m ] - E[ y | b_{l-1} < x_k \leq b_l ] \nonumber \\ & \nonumber \\ \approx & \frac{ \exp{( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_j x_j + \beta_k E[ x_k | b_{m-1} < x_k \leq b_m ])}} { 1 + \exp{( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_j x_j + \beta_k E[ x_k | b_{m-1} < x_k \leq b_m ] )}} \nonumber \\ & \nonumber \\ & -\frac{ \exp{( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_j x_j + \beta_k E[ x_k | b_{l-1} < x_k \leq b_l ])}} { 1 + \exp{( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_j x_j + \beta_k E[ x_k | b_{l-1} < x_k \leq b_l ] )}} \label{eq:logEffect} \end{align} where again $E[ x_k | b_{m-1} < x_k \leq b_m ] = \bar{x}_{km}$ can be approximated by the mid-points of the intervals \begin{equation} \bar{x}_{km} \approx \frac{ b_{m-1} + b_m }{2} \; \forall \; m = 1, \ldots , M. \end{equation} For model specifications that additionally include a quadratic term of the explanatory variable, \begin{align} E_{k,ml} = & E[ y | b_{m-1} < x_k \leq b_m ] - E[ y | b_{l-1} < x_k \leq b_l ] \nonumber \\ & \nonumber \\ \approx & \frac{ \exp{( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus (k, k+1)} \beta_j x_j + \beta_k E[ x_k | b_{m-1} < x_k \leq b_m ] + \beta_{k+1} E[ x_k^2 | b_{m-1} < x_k \leq b_m ])}} { 1 + \exp{( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus (k, k+1) } \beta_j x_j + \beta_k E[ x_k | b_{m-1} < x_k \leq b_m ] + \beta_{k+1} E[ x_k^2 | b_{m-1} < x_k \leq b_m ])}} \nonumber \\ & \nonumber \\ & -\frac{ \exp{( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus (k, k+1) } \beta_j x_j + \beta_k E[ x_k | b_{l-1} < x_k \leq b_l ] + \beta_{k+1} E[ x_k^2 | b_{l-1} < x_k \leq b_l ])}} { 1 + \exp{( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus (k, k+1) } \beta_j x_j + \beta_k E[ x_k | b_{l-1} < x_k \leq b_l ] + \beta_{k+1} E[ x_k^2 | b_{l-1} < x_k \leq b_l ])}} \label{eq:logQuadEffect} \end{align} where the values for $E[ x_k | b_{m-1} < x_k \leq b_m ]$ and $E[ x_k^2 | b_{m-1} < x_k \leq b_m ]$ remain the same as outlined in section~\ref{sec:linEffInt}. In the case of the multi-nomial logit function equation~(\ref{eq:logEffect}) modifies to \begin{align} E_{k,ml,p} \approx & \dfrac{ \exp{( \beta_{0,p} + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_{j,p} x_j + \beta_{k,p} E[ x_k | b_{m-1} < x_k \leq b_m ])}} { 1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp{( \beta_{0,p} + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_{j,p} x_j + \beta_{k,p} E[ x_k | b_{m-1} < x_k \leq b_m ] )}} \nonumber \\ & \nonumber \\ & -\dfrac{ \exp{( \beta_{0,p} + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_{j,p} x_j + \beta_{k,p} E[ x_k | b_{l-1} < x_k \leq b_l ])}} { 1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp{( \beta_{0,p} + \sum_{j \in \{ 1, \ldots, K \} \setminus k } \beta_{j,p} x_j + \beta_{k,p} E[ x_k | b_{l-1} < x_k \leq b_l ] )}} \label{eq:MNlogEffect} \end{align} and equation~(\ref{eq:logQuadEffect}) modifies accordingly. In the case of the conditional logit function equation~(\ref{eq:logEffect}) modifies to \begin{align} E_{k,ml,p} \approx & \dfrac{ \exp{( \gamma_0 + \sum_{j \in \{ 1, \ldots, T \} \setminus t } \gamma_j z_{j,p} + \gamma_k E[ z_{k,p} | b_{m-1} < z_{k,p} \leq b_m ])}} { \sum_{p \in C } \exp{( \gamma_0 + \sum_{j \in \{ 1, \ldots, T \} \setminus t } \gamma_j z_{j,p} + \gamma_k E[ z_{k,p} | b_{m-1} < z_{k,p} \leq b_m ] )}} \nonumber \\ & \nonumber \\ & -\dfrac{ \exp{( \gamma_0 + \sum_{j \in \{ 1, \ldots, T \} \setminus t } \gamma_j z_{j,p} + \gamma_k E[ z_{k,p} | b_{l-1} < z_{k,p} \leq b_l ])}} { \sum_{p \in C } \exp{( \gamma_0 + \sum_{j \in \{ 1, \ldots, T \} \setminus t } \gamma_j z_{j,p} + \gamma_k E[ z_{k,p} | b_{l-1} < z_{k,p} \leq b_l ] )}} \label{eq:CondlogEffect} \end{align} and equation~(\ref{eq:logQuadEffect}) modifies accordingly. And in the case of the nested logit function equation~(\ref{eq:logEffect}) modifies to \begin{align} E_{k, ml, p, o} =& \; \frac{ \exp( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \beta_k/\lambda_o E[ x_{k,p} | b_{m-1} < x_{k,p} \leq b_m ] ) } { \sum_{p \in B_o} \exp( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} ) + \beta_k/\lambda_o E[ x_{k,p} | b_{m-1} < x_{k,p} \leq b_m ] ) } \nonumber \\ & \cdot \frac{ \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_o )} { \sum_{o = 1}^O \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_o )} \nonumber \\ &- \frac{ \exp( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \beta_k/\lambda_o E[ x_{k,p} | b_{l-1} < x_{k,p} \leq b_l ] ) } { \sum_{p \in B_o} \exp( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} ) + \beta_k/\lambda_o E[ x_{k,p} | b_{l-1} < x_{k,p} \leq b_l ] ) } \nonumber \\ & \cdot \frac{ \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_o )} { \sum_{o = 1}^O \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_o )} \end{align} where equation~(\ref{eq:IVo}) modifies to \begin{equation} IV_o = \ln \sum_{p \in B_o} \exp \left( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \beta_k/\lambda_0 E[ x_{k,p}| b_{m-1} < x_{k,p} \leq b_m ] \right) \end{equation} and \begin{equation} IV_o = \ln \sum_{p \in B_o} \exp \left( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \beta_k/\lambda_0 E[ x_{k,p}| b_{l-1} < x_{k,p} \leq b_l ] \right), \end{equation} respectively. In order to calculate the standard error of $E_{k,ml}$, we again apply a simplified Delta method \begin{equation} \se( E_{k,lm}) = \sqrt{ \left( \frac{\partial E_{k,lm}}{\partial \beta}\right)^\top \cdot \var( \beta ) \cdot \frac{\partial E_{k, lm}}{\partial \beta }} \label{eq:logDelta3.3} \end{equation} where the elements of the gradient vector, $\dfrac{\partial E_{k, ml}}{\partial \mathbf{\beta}}$, are: \begin{align} \frac{\partial E_{k, ml}}{\partial \beta_0} =& \dfrac{ \exp_m( \cdot )}{[ 1 + \exp_m( \cdot )]^2 } - \dfrac{ \exp_l( \cdot )}{[ 1 + \exp_l( \cdot )]^2 } \label{eq:Sgrad3.3a} \\[1em] \frac{\partial E_{k, ml}}{\partial \beta_j } =& \dfrac{ \exp_m( \cdot ) \cdot \bar{x}_j }{[ 1 + \exp_m( \cdot )]^2 } - \dfrac{ \exp_l( \cdot ) \cdot \bar{x}_j }{[ 1 + \exp_l( \cdot )]^2 } \;\;\; \forall \;\; \beta_j \;\; j \neq k, k + 1 \\[1em] \frac{\partial E_{k, ml}}{\partial \beta_k} =& \dfrac{ \exp_m( \cdot ) \cdot \bar{x}_{km}}{[ 1 + \exp_m( \cdot )]^2 } - \dfrac{ \exp_l( \cdot ) \cdot \bar{x}_{kl}}{[ 1 + \exp_l( \cdot )]^2 } \\[1em] \frac{\partial E_{k, ml}}{\partial \beta_{k+1}} =& \dfrac{ \exp_m( \cdot ) \cdot \overline{x^2}_{km}}{[ 1 + \exp_m( \cdot )]^2 } - \dfrac{ \exp_l( \cdot ) \cdot \overline{x^2}_{kl}}{[ 1 + \exp_l( \cdot )]^2 } \label{eq:Sgrad3.3b} \end{align} where \begin{equation} \exp_m(\cdot) \equiv \exp{\left( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus (k, k+1)} \beta_j x_j + \beta_k E[ x_k | b_{m-1} < x_k \leq b_m ] + \beta_{k+1} E[ x_k^2 | b_{m-1} < x_k \leq b_m ]\right)} \nonumber \end{equation} and \begin{equation} \exp_l(\cdot) \equiv \exp{\left( \beta_0 + \sum_{j \in \{ 1, \ldots, K \} \setminus (k, k+1)} \beta_j x_j + \beta_k E[ x_k | b_{l-1} < x_k \leq b_l ] + \beta_{k+1} E[ x_k^2 | b_{l-1} < x_k \leq b_l ]\right)}. \nonumber \end{equation} In the case of the multi-nomial logit the gradient vector, $\dfrac{\partial E_{k, ml, p}}{\partial \mathbf{\beta}}$, modifies to: \begin{align} \frac{ \partial E_{k, ml,p}}{\partial \beta_{0,p}} =& \, \frac{ \exp_{m, p}( \cdot ) \cdot ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ m, o } ( \cdot ) ) } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \label{eq:MNLgrad3.3a} \\ & - \frac{ \exp_{l, p}( \cdot ) \cdot ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ l, o } ( \cdot )) } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ l, p }( \cdot ) ]^2 } \nonumber \\[1em] \frac{ \partial E_{k, ml,p}}{\partial \beta_{0,o}} =& \, \frac{ \exp_{ l, p }( \cdot ) \exp_{ l, o }( \cdot )} {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ l, p }( \cdot ) ]^2 } \\ & - \frac{ \exp_{ m, p }( \cdot ) \exp_{ m, o }( \cdot )} {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \;\;\; \forall \; o \neq p \nonumber \\[1em] \frac{ \partial E_{k, ml,p}}{\partial \beta_{j,p}} =& \, \frac{ \exp_{ m, p }( \cdot ) ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ m, o } ( \cdot ) ) \cdot x_j } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \\ & - \frac{ \exp_{ l, p }( \cdot ) ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ l, o } ( \cdot ) ) \cdot x_j } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ l, p }( \cdot ) ]^2 } \;\;\; \forall \; j \in \{ 1, \ldots, K \} \setminus k, k + 1 \nonumber \\[1em] \frac{ \partial E_{k, ml,p}}{\partial \beta_{j,o}} =& \, \frac{ \exp_{ l, p }( \cdot ) \exp_{ l, o }( \cdot ) \cdot x_j } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ l, p }( \cdot ) ]^2 } \\ & - \frac{ \exp_{ m, p }( \cdot ) \exp_{ m, o }( \cdot ) \cdot x_j } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \;\;\; \forall \; o \neq p, \; j \in \{ 1, \ldots, K \} \setminus k, k+1 \nonumber \\[1em] \frac{ \partial E_{k, ml,p}}{\partial \beta_{k,p}} =& \, \frac{ \exp_{ m, p }( \cdot ) ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ m, o } ( \cdot ) ) \cdot \bar{x}_k } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \\ & - \frac{ \exp_{ l, p }( \cdot ) ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ l, o } ( \cdot ) ) \cdot \bar{x}_k } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ l, p }( \cdot ) ]^2 } \nonumber \\[1em] \frac{ \partial E_{k, ml,p}}{\partial \beta_{k,o}} =& \, \frac{ \exp_{l,p}( \cdot ) \exp_{ l, o }( \cdot ) \cdot \bar{x}_k } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ l, p }( \cdot ) ]^2 } - \frac{ \exp_{ m, p }( \cdot ) \exp_{ m, o }( \cdot ) \cdot \bar{x}_k } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \\[1em] \frac{ \partial E_{k, ml,p}}{\partial \beta_{k+1,p}} =& \, \frac{ \exp_{ m, p }( \cdot ) ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ m, o } ( \cdot ) ) \cdot \overline{x^2_k} } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \\ & - \frac{ \exp_{ l, p }( \cdot ) ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ l, o } ( \cdot ) ) \cdot \overline{x^2_k} } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ l, p }( \cdot ) ]^2 } \nonumber \\[1em] \frac{ \partial E_{k, ml,p}}{\partial \beta_{k+1,o}} =& \, \frac{ \exp_{ l, p }( \cdot ) \exp_{ l, o }( \cdot ) \cdot \overline{x^2}_{km} } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ l, p }( \cdot ) ]^2 } - \frac{ \exp_{ m, p }( \cdot ) \exp_{ m, o }( \cdot ) \cdot \overline{x^2}_{km} } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ m, p }( \cdot ) ]^2 } \label{eq:MNLgrad3.3b} \end{align} where \begin{equation} \exp_{m, p}(\cdot) \equiv \exp{\left( \beta_{0,p} + \sum_{j \in \{ 1, \ldots, K \} \setminus (k, k+1)} \beta_{j,p} x_j + \beta_{k,p} E[ x_k | b_{m-1} < x_k \leq b_m ] + \beta_{k+1,p} E[ x_k^2 | b_{m-1} < x_k \leq b_m ]\right)} \nonumber \end{equation} and \begin{equation} \exp_{l, p}(\cdot) \equiv \exp{\left( \beta_{0,p} + \sum_{j \in \{ 1, \ldots, K \} \setminus (k, k+1)} \beta_{j,p} x_j + \beta_{k,p} E[ x_k | b_{l-1} < x_k \leq b_l ] + \beta_{k+1,p} E[ x_k^2 | b_{l-1} < x_k \leq b_l ]\right)}. \nonumber \end{equation} In the case of the conditional logit the gradient vector, $\dfrac{\partial E_{k, ml, p}}{\partial \mathbf{\beta}}$, modifies to: \begin{align} \frac{\partial E_{t, ml, p}}{\partial \gamma_0} =& \frac{ \exp_{m,p}(\cdot) \cdot \sum_{p\in C} \exp_{m,p}(\cdot) - \exp_{m,p}(\cdot) \cdot \sum_{p\in C} \exp_{m,p}(\cdot) } { \left( \sum_{p\in C} \exp_{m,p}(\cdot) \right)^2 } \label{eq:condLoggrad3.3a}\\ & - \frac{ \exp_{l,p}(\cdot) \cdot \sum_{p\in C} \exp_{l,p}(\cdot) - \exp_{l,p}(\cdot) \cdot \sum_{p\in C} \exp_{l,p}(\cdot) } { \left( \sum_{p\in C} \exp_{l,p}(\cdot) \right)^2 } \nonumber \\[1em] \frac{\partial E_{t, ml, p}}{\partial \gamma_j } =& \frac{ \exp_{m,p}(\cdot) z_{j,p} \cdot \sum_{p\in C} \exp_{m,p}(\cdot) - \exp_{m,p}(\cdot) \cdot \sum_{p\in C} (\exp_{m,p}(\cdot) z_{j,p} ) } { \left( \sum_{p\in C} \exp_{m,p}(\cdot) \right)^2 } \\ & - \frac{ \exp_{l,p}(\cdot) z_{j,p} \cdot \sum_{p\in C} \exp_{l,p}(\cdot) - \exp_{l,p}(\cdot) \cdot \sum_{p\in C} ( \exp_{l,p}(\cdot) z_{j,p} ) } { \left( \sum_{p\in C} \exp_{l,p}(\cdot) \right)^2 } \;\;\; \forall \;\; \gamma_j \;\; j \neq t, t + 1 \nonumber \\[1em] \frac{\partial E_{t, ml, p}}{\partial \gamma_t} =& \frac{ \exp_{m,p}(\cdot) z_{t,p} \cdot \sum_{p\in C} \exp_{m,p}(\cdot) - \exp_{m,p}(\cdot) \cdot \sum_{p\in C} ( \exp_{m,p}(\cdot) z_{t,p} ) } { \left( \sum_{p\in C} \exp_{m,p}(\cdot) \right)^2 } \\ & - \frac{ \exp_{l,p}(\cdot) z_{t,p} \cdot \sum_{p\in C} \exp_{l,p}(\cdot) - \exp_{l,p}(\cdot) \cdot \sum_{p\in C} ( \exp_{l,p}(\cdot) z_{t,p} ) } { \left( \sum_{p\in C} \exp_{l,p}(\cdot) \right)^2 } \nonumber \\[1em] \frac{\partial E_{t, ml, p}}{\partial \gamma_{t+1}} =& \frac{ \exp_{m,p}(\cdot) z_{t+1,p} \cdot \sum_{p\in C} \exp_{m,p}(\cdot) - \exp_{m,p}(\cdot) \cdot \sum_{p\in C} ( \exp_{m,p}(\cdot) z_{t+1,p} ) } { \left( \sum_{p\in C} \exp_{m,p}(\cdot) \right)^2 } \label{eq:condLoggrad3.3b} \\ & - \frac{ \exp_{l,p}(\cdot) z_{t+1,p} \cdot \sum_{p\in C} \exp_{l,p}(\cdot) - \exp_{l,p}(\cdot) \cdot \sum_{p\in C} ( \exp_{l,p}(\cdot) z_{t+1,p} ) } { \left( \sum_{p\in C} \exp_{l,p}(\cdot) \right)^2 } \nonumber \end{align} where \begin{equation} \exp_{m, p}(\cdot) \equiv \exp{\left( \gamma_0 + \sum_{j \in \{ 1, \ldots, T \} \setminus (t, t+1)} \gamma_j z_{j,p} + \gamma_{t,p} E[ z_t | b_{m-1} < z_t \leq b_m ] + \gamma_{t+1,p} E[ z_t^2 | b_{m-1} < z_t \leq b_m ]\right)} \nonumber \end{equation} and \begin{equation} \exp_{l, p}(\cdot) \equiv \exp{\left( \gamma_0 + \sum_{j \in \{ 1, \ldots, T \} \setminus (t, t+1)} \gamma_j z_{j,p} + \gamma_t E[ z_t | b_{l-1} < z_t \leq b_l ] + \gamma_t+1 E[ z_t^2 | b_{l-1} < z_t \leq b_l ]\right)}. \nonumber \end{equation} Finally, in the case of the nested logit, the gradient vector, $\dfrac{\partial E_{k, ml, p}}{\partial \mathbf{\beta}}$, modifies to: \begin{align} \frac{\partial E_{t, ml, p}}{\partial \beta_0} =& \left( \left( \frac{1}{\lambda_o} - \frac{\sum_{p \in B_o} \left( \exp_{m,p,o} (\cdot) \frac{1}{\lambda_o} \right)} {\sum_{p \in B_o} \exp_{m,p,o} (\cdot) } \right) \right. \nonumber \\ & \left. + \left( \sum_{p \in B_o} 1 - \frac{ \sum_{o=1}^O \left( \exp_{m,o} (\cdot) \sum_{p \in B_o} 1 \right)} { \sum_{o=1}^O \exp_{m,o} (\cdot)} \right) \right) \pi_{m,p} \pi_{m,o} \nonumber \\ & - \left( \left( \frac{1}{\lambda_o} - \frac{\sum_{p \in B_o} \left( \exp_{l,p,o} (\cdot) \frac{1}{\lambda_o} \right)} {\sum_{p \in B_o} \exp_{l,p,o} (\cdot) } \right) \right. \nonumber \\ & \left. + \left( \sum_{p \in B_o} 1 - \frac{ \sum_{o=1}^O \left( \exp_{l,o} (\cdot) \sum_{p \in B_o} 1 \right)} { \sum_{o=1}^O \exp_{l,o} (\cdot)} \right) \right) \pi_{l,p} \pi_{l,o} \\[1em] %========================================================================================= \frac{\partial E_{t, ml, p}}{\partial \beta_j} =& \left( \left( \frac{x_{p,j}}{\lambda_o} - \frac{\sum_{p \in B_o} \left( \exp_{m,p,o} (\cdot) \frac{x_{p,j}}{\lambda_o} \right)} {\sum_{p \in B_o} \exp_{m,p,o} (\cdot) } \right) \right. \nonumber \\ & \left. + \left( \sum_{p \in B_o} x_{p,j} - \frac{ \sum_{o=1}^O \left( \exp_{m,o} (\cdot) \sum_{p \in B_o} x_{p,j} \right)} { \sum_{o=1}^O \exp_{m,o} (\cdot)} \right) \right) \pi_{m,p} \pi_{m,o} \nonumber \\ & - \left( \left( \frac{x_{p,j}}{\lambda_o} - \frac{\sum_{p \in B_o} \left( \exp_{l,p,o} (\cdot) \frac{x_{p,j}}{\lambda_o} \right)} {\sum_{p \in B_o} \exp_{l,p,o} (\cdot) } \right) \right. \nonumber \\ & \left. + \left( \sum_{p \in B_o} x_{p,j} - \frac{ \sum_{o=1}^O \left( \exp_{l,o} (\cdot) \sum_{p \in B_o} x_{p,j} \right)} { \sum_{o=1}^O \exp_{l,o} (\cdot)} \right) \right) \pi_{l,p} \pi_{l,o} \\[1em] %========================================================================================= \frac{\partial E_{t, ml, p}}{\partial \beta_k} =& \left( \left( \frac{x_{p,m}}{\lambda_o} - \frac{\sum_{p \in B_o} \left( \exp_{m,p,o} (\cdot) \frac{x_{p,m}}{\lambda_o} \right)} {\sum_{p \in B_o} \exp_{m,p,o} (\cdot) } \right) \right. \nonumber \\ & \left. + \left( \sum_{p \in B_o} x_{p,m} - \frac{ \sum_{o=1}^O \left( \exp_{m,o} (\cdot) \sum_{p \in B_o} x_{p,m} \right)} { \sum_{o=1}^O \exp_{m,o} (\cdot)} \right) \right) \pi_{m,p} \pi_{m,o} \nonumber \\ & - \left( \left( \frac{x_{p,l}}{\lambda_o} - \frac{\sum_{p \in B_o} \left( \exp_{l,p,o} (\cdot) \frac{x_{p,l}}{\lambda_o} \right)} {\sum_{p \in B_o} \exp_{l,p,o} (\cdot) } \right) \right. \nonumber \\ & \left. + \left( \sum_{p \in B_o} x_{p,l} - \frac{ \sum_{o=1}^O \left( \exp_{l,o} (\cdot) \sum_{p \in B_o} x_{p,l} \right)} { \sum_{o=1}^O \exp_{l,o} (\cdot)} \right) \right) \pi_{l,p} \pi_{l,o} \end{align} \begin{align} \frac{\partial E_{t, ml, p}}{\partial \beta_{k+1}} =& \left( \left( \frac{x_{p,m+1}}{\lambda_o} - \frac{\sum_{p \in B_o} \left( \exp_{m,p,o} (\cdot) \frac{x_{p,m+1}}{\lambda_o} \right)} {\sum_{p \in B_o} \exp_{m,p,o} (\cdot) } \right) \right. \nonumber \\ & \left. + \left( \sum_{p \in B_o} x_{p,m+1} - \frac{ \sum_{o=1}^O \left( \exp_{m,o} (\cdot) \sum_{p \in B_o} x_{p,m+1} \right)} { \sum_{o=1}^O \exp_{m,o} (\cdot)} \right) \right) \pi_{m,p} \pi_{m,o} \nonumber \\ & - \left( \left( \frac{x_{p,l+1}}{\lambda_o} - \frac{\sum_{p \in B_o} \left( \exp_{l,p,o} (\cdot) \frac{x_{p,l+1}}{\lambda_o} \right)} {\sum_{p \in B_o} \exp_{l,p,o} (\cdot) } \right) \right. \nonumber \\ & \left. + \left( \sum_{p \in B_o} x_{p,l+1} - \frac{ \sum_{o=1}^O \left( \exp_{l,o} (\cdot) \sum_{p \in B_o} x_{p,l+1} \right)} { \sum_{o=1}^O \exp_{l,o} (\cdot)} \right) \right) \pi_{l,p} \pi_{l,o} \end{align} with \begin{align} \exp_{m,p,o}(\cdot) &\equiv \exp \left( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \beta_k/\lambda_o E[ x_{k,p} | b_{m-1} < x_{k,p} \leq b_m ] \right) \\[1em] \exp_{l,p,o}(\cdot) &\equiv \exp \left( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \beta_k/\lambda_o E[ x_{k,p} | b_{m-1} < x_{k,p} \leq b_m ] \right) \\[1em] \exp_{m/l,o} (\cdot) &\equiv \exp \left( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_o \right) \end{align} Using helper functions \code{EXSquared} (equation~\ref{eq:linEffectXSquaredBar}), \code{elaIntBounds} (checking arguments \code{refBound} and \code{intBound}), \code{checkXPos} (checking argument \code{xPos}, see appendix~\ref{sec:helperFunctions}), and \code{checkXBeta} (checking $\beta_0 + \sum_{j=1}^K \beta_j x_j$ for plausible values, see appendix~\ref{sec:helperFunctions}), the following function calculates the effect and its standard error according to equations~(\ref{eq:logEffect}), (\ref{eq:logQuadEffect}), and (\ref{eq:Sgrad3.3a}) to (\ref{eq:Sgrad3.3b}) for the binary logit, according to equations~(\ref{eq:MNlogEffect}) and (\ref{eq:MNLgrad3.3a}) to (\ref{eq:MNLgrad3.3b}) for the multi-nomial logit, and according to equations~(\ref{eq:CondlogEffect}) and (\ref{eq:condLoggrad3.3a}) to (\ref{eq:condLoggrad3.3b}) for the conditional logit: << >>= logEffInt <- function( allCoef, allCoefBra = NA, allXVal, allXValBra=NA, xPos, refBound, intBound, yCat, yCatBra, lambda, allCoefSE = rep( NA, length( allCoef ) ), method = "binary" ){ if( method == "binary" ){ # number of coefficients nCoef <- length( allCoef ) # check arguments if( length( allXVal ) != nCoef ){ stop( "argument 'allCoef' and 'allXVal' must have the same length" ) } if( length( allCoefSE ) != nCoef ){ stop( "argument 'allCoef' and 'allCoefSE' must have the same length" ) } } else if( method == "MNL" ){ # number of coefficients NCoef <- length( allCoef ) mCoef <- matrix( allCoef, nrow = length( allXVal )) nCoef <- dim( mCoef )[1] pCoef <- dim( mCoef )[2] # check arguments if( length( allXVal ) != nCoef ){ stop( "argument 'allCoef' and 'allXVal' must have the same length" ) } if( length( allCoefSE ) != NCoef ){ stop( "argument 'allCoef' and 'allCoefSE' must have the same length" ) } } else if( method == "CondL"){ # number of coefficients nCoef <- length( allCoef ) mXVal <- matrix( allXVal, nrow = nCoef ) pCoef <- dim( mXVal )[2] # check arguments if( dim( mXVal )[1] != nCoef ){ stop( "argument 'allCoef' and 'allXVal' must have the same length" ) } if( length( allCoefSE ) != nCoef ){ stop( "argument 'allCoef' and 'allCoefSE' must have the same length" ) } } else{ nCoef <- length( allCoef ) NCoef <- length( allCoefBra ) mXValBra <- matrix( allXValBra, nrow = NCoef ) nXValBra <- dim( mXValBra )[1] pXValBra <- dim( mXValBra )[2] # check arguments if( NCoef != nXValBra ){ stop( "arguments 'allCoefBra' and 'allXValBra' must have the same length") } O <- length( allXVal ) nXVal <- unlist( lapply( allXVal, function(x) dim( x )[1] ) ) pCoef <- unlist( lapply( allXVal, function(x) dim( x )[2] ) ) if( nCoef != nXVal[ yCatBra ] ){ stop( "arguments 'allCoef' and 'allXVal' must have the same length" ) } if( nCoef != length( allCoefSE) ){ stop( "arguments 'allCoef' and 'allCoefSE' must have the same length" ) } } checkXPos( xPos, minLength = 1, maxLength = 2, minVal = 1, maxVal = nCoef ) refBound <- elaIntBounds( refBound, 1, argName = "refBound" ) intBound <- elaIntBounds( intBound, 1, argName = "intBound" ) if( method == "binary" || method == "MNL" ){ if( any( !is.na( allXVal[ xPos ] ) ) ) { allXVal[ xPos ] <- NA warning( "values of argument 'allXVal[ xPos ]' are ignored", " (set these values to 'NA' to avoid this warning)" ) } }else if( method == "CondL" ){ for( p in 1:pCoef ){ if( any( !is.na( mXVal[ xPos, p ] ) ) ){ mXVal[ xPos, p ] <- NA warning( "values of argument 'allXVal[ xPos ]' are ignored", " (set these values to 'NA' to avoid this warning)" ) } } }else{ for( p in 1:pCoef[ yCatBra ] ){ if( any( !is.na( allXVal[[ yCatBra ]][ xPos, p ] ) ) ){ mXVal[ xPos, p ] <- NA warning( "values of argument 'allXVal[ xPos ]' are ignored", " (set these values to 'NA' to avoid this warning)" ) } } } # calculate xBars intX <- mean( intBound ) refX <- mean( refBound ) if( length( xPos ) == 2 ) { intX <- c( intX, EXSquared( intBound[1], intBound[2] ) ) refX <- c( refX, EXSquared( refBound[1], refBound[2] ) ) } if( length( intX ) != length( xPos ) || length( refX ) != length( xPos ) ) { stop( "internal error: 'intX' or 'refX' does not have the expected length" ) } # define X' * beta if( method == "binary" ){ intXbeta <- sum( allCoef * replace( allXVal, xPos, intX ) ) refXbeta <- sum( allCoef * replace( allXVal, xPos, refX ) ) checkXBeta( c( intXbeta, refXbeta ) ) } else if( method == "MNL" ){ intXbeta <- colSums( mCoef * replace( allXVal, xPos, intX ) ) refXbeta <- colSums( mCoef * replace( allXVal, xPos, refX ) ) } else if( method == "CondL" ){ mXValint <- mXValref <- mXVal for( p in 1:pCoef ){ mXValint[ ,p] <- replace( mXValint[ ,p], xPos, intX ) mXValref[ ,p] <- replace( mXValref[ ,p], xPos, refX ) } intXbeta <- colSums( allCoef * mXValint ) refXbeta <- colSums( allCoef * mXValref ) } else{ mCoef <- matrix( rep( allCoef, O ), nrow = nCoef, O ) %*% diag( 1/ lambda ) mXValint <- mXValref <- allXVal for( i in 1:O ){ for( p in 1:pCoef[i] ){ mXValint[[i]][ ,p] <- replace( mXValint[[i]][ ,p], xPos, intX ) mXValref[[i]][ ,p] <- replace( mXValref[[i]][ ,p], xPos, refX ) } } refXbeta <- intXbeta <- rep( list( NA ), O ) for( l in 1:O ){ intXbeta[[ l ]] <- colSums( mCoef[ ,l ] * mXValint[[ l ]] ) refXbeta[[ l ]] <- colSums( mCoef[ ,l ] * mXValref[[ l ]] ) } XbetaBra <- colSums( allCoefBra * mXValBra ) } # effect E_{k,ml} if( method == "binary" ){ eff <- exp( intXbeta )/( 1 + exp( intXbeta ) ) - exp( refXbeta )/( 1 + exp( refXbeta ) ) } else if( method == "MNL" ){ eff <- exp( intXbeta[ yCat ] )/( 1 + sum( exp( intXbeta ) ) ) - exp( refXbeta[ yCat ] )/( 1 + sum( exp( refXbeta ) ) ) } else if( method == "CondL"){ eff <- exp( intXbeta[ yCat ] )/( sum( exp( intXbeta ) ) ) - exp( refXbeta[ yCat ] )/( sum( exp( refXbeta ) ) ) } else{ intBranch <- refBranch <- rep( list( NA ), O ) for( l in 1:O ){ intBranch[[ l ]] <- exp( XbetaBra[ l ] + lambda[ l ] * log( sum( exp( intXbeta[[ l ]] ) ) ) ) refBranch[[ l ]] <- exp( XbetaBra[ l ] + lambda[ l ] * log( sum( exp( refXbeta[[ l ]] ) ) ) ) } intBranch <- unlist( intBranch ) refBranch <- unlist( refBranch ) eff <- exp( intXbeta[[ yCatBra ]][ yCat ] )/( sum( exp( intXbeta[[ yCatBra ]] ) ) ) * intBranch[ yCatBra ]/ sum( intBranch ) - exp( refXbeta[[ yCatBra ]][ yCat ] )/( sum( exp( refXbeta[[ yCatBra ]] ) ) ) * refBranch[ yCatBra ]/ sum( refBranch ) } # calculating approximate standard error # partial derivative of E_{k,ml} w.r.t. all estimated coefficients if( method == "binary" ){ derivCoef <- rep( NA, nCoef ) derivCoef[ -xPos ] <- ( exp( intXbeta )/( 1 + exp( intXbeta ) )^2 - exp( refXbeta )/( 1 + exp( refXbeta ) )^2 ) * allXVal[ -xPos ] derivCoef[ xPos ] <- exp( intXbeta )/( 1 + exp( intXbeta ) )^2 * intX - exp( refXbeta )/( 1 + exp( refXbeta ) )^2 * refX } else if( method == "MNL" ){ derivCoef <- matrix( NA, nrow=nCoef, ncol=pCoef ) for( p in 1:pCoef ){ if( p == yCat ){ derivCoef[ -xPos, p ] <- ( exp( intXbeta[ p ] ) * ( 1 + sum( exp( intXbeta[ -yCat ] ) ) )/ ( 1 + sum( exp( intXbeta ) ) )^2 - exp( refXbeta[ p ] ) * ( 1 + sum( exp( refXbeta[ -yCat ] ) ) )/ ( 1 + sum( exp( refXbeta ) ) )^2 ) * allXVal[ - xPos ] derivCoef[ xPos, p ] <- ( exp( intXbeta[ p ] ) * ( 1 + sum( exp( intXbeta[ -yCat ] ) ) )/ ( 1 + sum( exp( intXbeta ) ) )^2 ) * intX - ( exp( refXbeta[ p ] ) * ( 1 + sum( exp( refXbeta[ -yCat ] ) ) )/ ( 1 + sum( exp( refXbeta ) ) )^2 ) * refX } else{ derivCoef[ -xPos, p ] <- ( ( exp( refXbeta[ yCat ] ) * exp( refXbeta[ p ] ) )/ ( 1 + sum( exp( refXbeta ) ) )^2 - ( exp( intXbeta[ yCat ] ) * exp( intXbeta[ p ] ) )/ ( 1 + sum( exp( intXbeta ) ) )^2 ) * allXVal[ -xPos ] derivCoef[ xPos, p ] <- ( ( exp( refXbeta[ yCat ] ) * exp( refXbeta[ p ] ) )/ ( 1 + sum( exp( refXbeta ) ) )^2 ) * intX - ( ( exp( intXbeta[ yCat ] ) * exp( intXbeta[ p ] ) )/ ( 1 + sum( exp( intXbeta ) ) )^2 ) * refX } } derivCoef <- c( derivCoef ) } else if( method == "CondL" ){ derivCoef <- rep( NA, nCoef ) derivCoef[ -xPos ] <- ( exp( intXbeta[ yCat] ) * mXVal[ -xPos, yCat] * sum( exp( intXbeta ) ) - exp( intXbeta[ yCat] ) * rowSums( exp( intXbeta ) * mXVal[ -xPos, ] ) )/ ( sum( exp( intXbeta ) ) )^2 - ( exp( refXbeta[ yCat] ) * mXVal[ -xPos, yCat] * sum( exp( refXbeta ) ) - exp( refXbeta[ yCat] ) * rowSums( exp( refXbeta ) * mXVal[ -xPos, ] ) )/ ( sum( exp( refXbeta ) ) )^2 derivCoef[ xPos ] <- ( exp( intXbeta[ yCat] ) * intX * sum( exp( intXbeta ) ) - exp( intXbeta[ yCat] ) * sum( exp( intXbeta ) * intX ) )/ ( sum( exp( intXbeta ) ) )^2 - ( exp( refXbeta[ yCat] ) * refX * sum( exp( refXbeta ) ) - exp( refXbeta[ yCat] ) * sum( exp( refXbeta ) * refX ) )/ ( sum( exp( refXbeta ) ) )^2 } else{ derivCoef <- rep( NA, nCoef ) PImp <- exp( intXbeta[[ yCatBra ]][ yCat ])/( sum( exp( intXbeta[[ yCatBra ]] ) ) ) PIlp <- exp( refXbeta[[ yCatBra ]][ yCat ])/( sum( exp( refXbeta[[ yCatBra ]] ) ) ) PImo <- intBranch[ yCatBra ]/ sum( intBranch ) PIlo <- refBranch[ yCatBra ]/ sum( refBranch ) Om <- matrix( unlist( lapply( allXVal, function(x) rowSums( x[ -xPos, , drop = FALSE ] ) ) ), ncol = O ) derivCoef[ -xPos ] <- ( ( allXVal[[ yCatBra ]][ -xPos, yCat ]/lambda[ yCatBra ] - ( rowSums( ( allXVal[[ yCatBra ]][ -xPos, ]/lambda[ yCatBra ] ) %*% diag( exp( intXbeta[[ yCatBra ]] ) ) ) )/ ( sum( exp( intXbeta[[ yCatBra ]] ) ) ) ) + ( rowSums( allXVal[[ yCatBra ]][ -xPos, ] ) - ( rowSums( Om %*% diag( exp( intBranch ) ) )/ ( sum( intBranch ) ) ) ) ) * PImp * PImo - ( ( allXVal[[ yCatBra ]][ -xPos, yCat ]/lambda[ yCatBra ] - ( rowSums( ( allXVal[[ yCatBra ]][ -xPos, ]/lambda[ yCatBra ] ) %*% diag( exp( refXbeta[[ yCatBra ]] ) ) ) )/ ( sum( exp( refXbeta[[ yCatBra ]] ) ) ) ) + ( rowSums( allXVal[[ yCatBra ]][ -xPos, ] ) - ( rowSums( Om %*% diag( exp( refBranch ) ) )/ ( sum( refBranch ) ) ) ) ) * PIlp * PIlo derivCoef[ xPos ] <- ( ( intX/lambda[ yCatBra ] - ( sum( intX/lambda[ yCatBra ] * exp( intXbeta[[ yCatBra ]] ) ) )/ ( sum( exp( intXbeta[[ yCatBra ]] ) ) ) ) + ( intX * pCoef[ yCatBra ] - ( sum( intX * exp( intBranch ) )/ ( sum( intBranch ) ) ) ) ) * PImp * PImo - ( ( refX/lambda[ yCatBra ] - ( sum( refX/lambda[ yCatBra ] * exp( refXbeta[[ yCatBra ]] ) ) )/ ( sum( exp( refXbeta[[ yCatBra ]] ) ) ) ) + ( refX * pCoef[ yCatBra ] - ( sum( refX * exp( refBranch ) )/ ( sum( refBranch ) ) ) ) ) * PImp * PImo } # variance covariance of the coefficients (covariances set to zero) vcovCoef <- diag( allCoefSE^2 ) # approximate standard error of the effect effSE <- drop( sqrt( t( derivCoef ) %*% vcovCoef %*% derivCoef ) ) # object to be returned result <- c( effect = eff, stdEr = effSE ) return( result ) } @ where \begin{itemize} \item arguments \code{allCoef} $= (\beta_0, \ldots, \beta_K )^\top$ a vector of all coefficients from a binary logit model, \code{allCoef} $= ( \gamma_0, \ldots, \gamma_T )^\top$ a vector of all coefficients from a conditional logit model, \code{allCoefBra} $= ( \beta_0, \ldots, \beta_K )^\top$, a vector of all coefficients from a nested logit model, where \code{allCoefBra}, the coefficients at the branch level, must always be included in combination with \code{allCoef}, the coefficients at the twig level, or \code{allCoef} $= ( \beta_{01}, \ldots, \beta_{K1}, \ldots, \beta_{0P}, \ldots, \beta_{KP} )$ a vector of the $P$ sets of coefficients from the multi-nomial logit regression which does \textbf{not} include any values for the reference category; \item argument \code{allXVal} $= ( 1, \ldots, \bar{x}_K )^\top$, a vector of all corresponding sample means and 1 for the intercept for the binary or multi-nomial logit, \code{allXVal} $= ( 1, \ldots, \bar{z}_{T,1}, \ldots, 1, \ldots, \bar{z}_{T,P})$ a vector of the $P$ sets of explanatory variables from the conditional logit, or \code{allXVal} $=((( 1, \ldots, \bar{x}_{K,p,o} ), \ldots, (1, \ldots, \bar{x}_{K,P,o}))^\top, \ldots, (( 1, \ldots, \bar{x}_{K,p,O}), \ldots, ( 1, \ldots, \bar{x}_{K,P,O}))^\top )^\top$, a list where the list elements are the corresponding matrices of the sample means for each nest at the twig level,% \footnote{The reason why we make an exception from the usual vector in the case of the nested logit model is that nests, $o$, can have different dimensions wrt $P$.} and which must always be combined with the corresponding values at the branch level in \code{allXVal}; \item argument \code{allCoefSE} $= (\se(\beta_0), \ldots, \se( \beta_K))^\top$, a vector of standard errors for all coefficients of a binary or conditional logit regression or of the standard errors of the $P$ sets of coefficients in a multi-nomial logit regression; \item argument \code{xPos} $= k+1$ or $(k+1, k+2)^\top$ a scalar or vector of two elements indicating the position of the coefficients of interest in vectors \code{allCoef} and \code{allXVal}; \item argument \code{refBound} $= ( b_{l-1}, b_l )$ indicates the boundaries of the reference interval; \item argument \code{intBound} $= ( b_{m-1}, b_m )$ indicates the boundaries of the interval of interest; \item argument \code{method = c("binary", "MNL", "CondL", "NestL" )} indicating the estimator used, where option \code{"CondL"} can also be used to calculate the semi-elasticity and standard error of a nested logit at the branch level and where option \code{"NestL"} estimates the semi-elasticity and standard error of the combined probabilities at the branch and twig level; \item argument \code{"lambda"} is a vector of the $O$ parameter of the correlation coefficient between the the branch specific alternatives of interest, it is only relevant under option \code{"NestL"} and should be set to $1$ in case the estimation coefficients in \code{allCoef} are already corrected by $\dfrac{ \beta_k }{ \lambda_o }$; \item finally, \code{yCat} a scalar identifying which of the $P$ output categories is of interest, only in combination with \code{method = "MNL"}, \code{"CondL"}, and \code{"NestL"} (for the twig level), and \code{yCatBra} at the branch level. \end{itemize} << >>= # Example eff6a <- logEffInt( allCoef = c( 0.33, 0.22, 0.05, 0.6 ), allXVal = c( 1, NA, 0.16, 0.13 ), xPos = 2, refBound = c( 8, 12 ), intBound = c( 13, 15 ) ) eff6a eff6b <- logEffInt( allCoef = c( 0.33, 0.22, 0.05, 0.6 ), allXVal = c( 1, NA, 0.16, 0.13 ), xPos = 2, refBound = c( 8, 12 ), intBound = c( 13, 15 ), allCoefSE = c( 0.003, 0.045, 0.007, 0.009 ) ) eff6b # Example eff7a <- logEffInt( allCoef = c( 0.33, 0.22, 0.05, 0.6 ), allXVal = c( 1, NA, NA, 0.0004 ), xPos = c( 2, 3 ), refBound = c( 8, 12 ), intBound = c( 13, 15 )) eff7a eff7b <- logEffInt( allCoef = c( 0.33, 0.22, 0.05, 0.6 ), allXVal = c( 1, NA, NA, 0.13 ), xPos = c( 2, 3 ), refBound = c( 8, 12 ), intBound = c( 13, 15 ), allCoefSE = c( 0.003, 0.045, 0.007, 0.009 ) ) eff7b #Example eff8a <- logEffInt( allCoef = c( 0.2, 0.3, 0.5, -0.2, 0.03, 0.6 ), allXVal = c( 1, NA, 0.12 ), xPos = 2, refBound = c( 8, 12 ), intBound = c( 13, 15 ), yCat = 2, method = "MNL" ) eff8a eff8b <- logEffInt( allCoef = c( 0.2, 0.3, 0.5, -0.2, 0.03, 0.6 ), allXVal = c( 1, NA, 0.12 ), xPos = 2, refBound = c( 8, 12 ), intBound = c( 13, 15 ), yCat = 2, allCoefSE = c( 0.003, 0.045, 0.007, 0.009, 0.0008, 0.9 ), method = "MNL" ) eff8b #Example eff9a <- logEffInt( allCoef = c( 0.2, 0.3, 0.5, -0.2, 0.03, 0.6 ), allXVal = c( 1, NA, NA ), xPos = c( 2, 3 ), refBound = c( 8, 12 ), intBound = c( 13, 15 ), yCat = 2, method = "MNL" ) eff9a eff9b <- logEffInt( allCoef = c( 0.2, 0.3, 0.5, -0.2, 0.03, 0.6 ), allXVal = c( 1, NA, NA ), xPos = c( 2, 3 ), refBound = c( 8, 12 ), intBound = c( 13, 15 ), yCat = 2, allCoefSE = c( 0.003, 0.045, 0.007, 0.009, 0.0008, 0.9 ), method = "MNL" ) eff9b #Example eff10a <- logEffInt( allCoef = c( 0.2, 0.3, 0.5, 0.091 ), allXVal = c( 1, NA, NA, 2.45, 1, NA, NA, 0.79 ), xPos = c( 2, 3 ), refBound = c( 8, 12 ), intBound = c( 13, 15 ), yCat = 2, method = "CondL" ) eff10a eff10b <- logEffInt( allCoef = c( 0.2, 0.3, 0.5, 0.091 ), allXVal = c( 1, NA, NA, 2.45, 1, NA, NA, 0.79 ), xPos = c( 2, 3 ), refBound = c( 8, 12 ), intBound = c( 13, 15 ), allCoefSE = c( 0.003, 0.045, 0.007, 0.009 ), yCat = 2, method = "CondL" ) eff10b # Example matrix1 <- matrix( c( 1, NA, 0.3, 0.09, 1, NA, 0.9, 1.8 ), nrow = 4 ) matrix2 <- matrix( c( 1, NA, 0.099, 0.211 ), nrow = 4 ) eff12a <- logEffInt( allCoefBra = c( 0.445, 0.03, -0.002 ), allCoef = c( 1.8, 0.005, -0.12, 0.8 ), allXValBra = c( 1, 3.3, 4.5, 1, 0.1, 0.987 ), allXVal = list( matrix1, matrix2 ), refBound = c( 0.5, 1.5 ), intBound = c( 2, 3.5 ), xPos = 2, yCatBra = 1, yCat = 2, lambda = c( 0.8, 1 ), method = "NestedL" ) eff12a matrix1 <- matrix( c( 1, NA, 0.3, 0.09, 1, NA, 0.9, 1.8 ), nrow = 4 ) matrix2 <- matrix( c( 1, NA, 0.099, 0.211 ), nrow = 4 ) eff12b <- logEffInt( allCoefBra = c( 0.445, 0.03, -0.002 ), allCoef = c( 1.8, 0.005, -0.12, 0.8 ), allXValBra = c( 1, 3.3, 4.5, 1, 0.1, 0.987 ), allXVal = list( matrix1, matrix2 ), allCoefSE = c( 0.003, 0.045, 0.007, 0.0032 ), refBound = c( 0.5, 1.5 ), intBound = c( 2, 3.5 ), xPos = 2, yCatBra = 1, yCat = 2, lambda = c( 0.8, 1 ), method = "NestedL" ) eff12b @ %===================================================================================== \subsection{Grouping and re-basing categorical variables} \label{sec:3.4} We consider a regression model \begin{align} \Pr( y = 1 | x ) &= \frac{ \exp \left( \beta_0 + \sum_{j\in\{1, \ldots, K \} \setminus k} \beta_j x_j + \sum_{m \in \{1,...,M \} \setminus m^*} \delta_m D_m \right)} { 1 + \exp \left( \beta_0 + \sum_{j\in\{1, \ldots, K \} \setminus k} \beta_j x_j + \sum_{m \in \{1,...,M \} \setminus m^*} \delta_m D_m \right)} \label{eq:logInterv}\\ D_m &= \begin{cases} 1 & \text{ if } x_k \in c_m \\ 0 & \text{ otherwise} \end{cases} \quad \forall \; m = 1, \ldots , M. \end{align} Like in section~\ref{sec:lingroup}, the $k$th xplanatory variable is coded into $M$ mutually exclusive categories $c_1, \ldots, c_M$, with $c_m \cap c_l = \emptyset \; \forall \; m \neq l$, and $D_1, \ldots, D_M$ the corresponding dummy variables. In the case of a binary logit regression, equation~(\ref{eq:linEffGr}) modifies to \begin{align} E_{k,lr} = \; & E[ y | x_k \in c_l^* ] - E[ y | x_k \in c_r^* ] \\ = \; & \frac{ \exp \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m E[ D_m | x_k \in c_l^* ] \right)} { 1 + \exp \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m E[ D_m | x_k \in c_l^* ] \right)} \\ & - \frac{ \exp \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^{M} \delta_m E[ D_m | x_k \in c_r^* ] \right)} { 1 + \exp \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^{M} \delta_m E[ D_m | x_k \in c_r^* ] \right)} \nonumber \\ = \; & \frac{ \exp \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m D_{ml} \right)} { 1 + \exp \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m D_{ml} \right)} \label{eq:logitEffGr} \\ & - \frac{ \exp \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^{M} \delta_m D_{mr} \right)} { 1 + \exp \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m D_{mr} \right)} \nonumber \end{align} where $D_{mn}$ is defined as in equation~(\ref{eq:linEffGrD}). Whilst for the case of the multi-nomial logit regression, equation~(\ref{eq:linEffGr}) modifies to \begin{align} E_{k, lr, p} = \;& \frac{ \exp \left( \beta_{0,p} + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_{j,p} x_j + \sum_{m=1}^M \delta_{m,p} D_{ml} \right)} { 1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp \left( \beta_{0,p} + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_{j,p} x_j + \sum_{m=1}^M \delta_{m,p} D_{ml} \right)} \label{eq:MNlogitEffGr} \\ & - \frac{ \exp \left( \beta_{0,p} + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_{j,p} x_j + \sum_{m=1}^{M} \delta_{m,p} D_{mr} \right)} { 1 + \sum_{p \in \{1, \ldots, P \} \setminus p^* } \exp \left( \beta_{0,p} + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_{j,p} x_j + \sum_{m=1}^M \delta_{m,p} D_{mr} \right)} \nonumber \end{align} and for the case of the conditional logit regression, equation~(\ref{eq:linEffGr}) modifies to \begin{align} E_{t, lr, p} = \;& \frac{ \exp \left( \gamma_0 + \sum_{j \in \{ 1, \ldots , T \} \setminus t} \gamma_j z_{j,p} + \sum_{m=1}^M \delta_m D_{ml} \right)} { \sum_{p \in C } \exp \left( \gamma_0 + \sum_{j \in \{ 1, \ldots , T \} \setminus t} \gamma_j z_{j,p} + \sum_{m=1}^M \delta_m D_{ml} \right)} \label{eq:CondlogitEffGr} \\ & - \frac{ \exp \left( \gamma_0 + \sum_{j \in \{ 1, \ldots , T \} \setminus t} \gamma_j z_{j,p} + \sum_{m=1}^{M} \delta_m D_{mr} \right)} { \sum_{p \in C } \exp \left( \gamma_0 + \sum_{j \in \{ 1, \ldots , T \} \setminus t} \gamma_j z_{j,p} + \sum_{m=1}^M \delta_m D_{mr} \right)} \nonumber \end{align} and for the case of the nested logit regression, equation~(\ref{eq:linEffGr}) modifies to \begin{align} E_{t, lr, p, o} = \;& \frac{ \exp( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \sum_{m=1}^M \delta_m/\lambda_o D_{ml} ) } { \sum_{p \in B_o} \exp( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \sum_{m=1}^M \delta_m/\lambda_o D_{ml} ) } \nonumber \\ & \cdot \frac{ \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_{o,l} ) } { \sum_{o = 1}^O \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_{o,l} ) } \nonumber \\ & - \frac{ \exp( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \sum_{m=1}^M \delta_m/\lambda_o D_{mr} ) } { \sum_{p \in B_o} \exp( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \sum_{m=1}^M \delta_m/\lambda_o D_{mr} ) } \nonumber \\ & \cdot \frac{ \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_{o,r} ) } { \sum_{o = 1}^O \exp( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_{o,r} ) } \end{align} where \begin{equation} IV_{o, l/r} = \ln \sum_{p \in B_o} \exp \left( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \sum_{m=1}^M \gamma_m/\lambda_0 D_{m, l/r} \right) . \end{equation} In order to calculate the approximate standard error of the new effect $E_{k,lr}$, we again apply the Delta method: \begin{equation} \se( E_{k,lr} ) = \sqrt{ \left( \frac{ \partial E_{k,lr}} { \partial \begin{pmatrix} \beta \\ \delta \end{pmatrix}} \right)^\top \cdot \var \begin{pmatrix} \beta \\ \delta \end{pmatrix} \cdot \frac{ \partial E_{k,lr}} { \partial \begin{pmatrix} \beta \\ \delta \end{pmatrix}}} , \label{eq:34Delta} \end{equation} with the partial derivatives defined as \begin{align} \frac{\partial E_{k,lr}}{\partial \beta_0} &= \frac{\exp_{ml}(\cdot)}{[1 + \exp_{ml}(\cdot)]^2} - \frac{\exp_{mr}(\cdot)}{[1 + \exp_{mr}(\cdot)]^2} \\ \frac{\partial E_{k,lr}}{\partial \beta_j} &= \left( \frac{\exp_{ml}(\cdot)}{[1 + \exp_{ml}(\cdot)]^2} - \frac{\exp_{mr}(\cdot)}{[1 + \exp_{mr}(\cdot)]^2}\right) \cdot \bar{x}_j \; \forall \; j \in \{ 1, \ldots , K \} \setminus k \\ \frac{\partial E_{k,lr}}{\partial \delta_m} &= \frac{\exp_{ml}(\cdot)}{[1 + \exp_{ml}(\cdot)]^2} D_{ml} - \frac{\exp_{mr}(\cdot)}{[1 + \exp_{mr}(\cdot)]^2} D_{mr} \; \forall \; m = 1, \ldots, M. \end{align} and for the multi-nomial logit as \begin{align} \frac{ \partial E_{k, lr,p}}{\partial \beta_{0,p}} =& \, \frac{ \exp_{ml, p}( \cdot ) \cdot ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ ml, o } ( \cdot ) ) } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ ml, p }( \cdot ) ]^2 } \label{eq:MNLGrad1} \\ & - \frac{ \exp_{mr, p}( \cdot ) \cdot ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ mr, o } ( \cdot )) } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ mr, p }( \cdot ) ]^2 } \nonumber \\[1em] \frac{ \partial E_{k, lr, p}}{\partial \beta_{0,o}} =& \, \frac{ \exp_{ mr, p }( \cdot ) \exp_{ mr, o }( \cdot )} {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ mr, p }( \cdot ) ]^2 } \\ & - \frac{ \exp_{ ml, p }( \cdot ) \exp_{ ml, o }( \cdot )} {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ ml, p }( \cdot ) ]^2 } \;\;\; \forall \; o \neq p \nonumber \\[1em] \frac{ \partial E_{k, lr,p}}{\partial \beta_{j,p}} =& \, \frac{ \exp_{ ml, p }( \cdot ) ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ ml, o } ( \cdot ) ) \cdot x_j } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ ml, p }( \cdot ) ]^2 } \\ & - \frac{ \exp_{ mr, p }( \cdot ) ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ mr, o } ( \cdot ) ) \cdot x_j } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ mr, p }( \cdot ) ]^2 } \;\;\; \forall \; j \in \{ 1, \ldots, K \} \setminus k, k + 1 \nonumber \\[1em] \frac{ \partial E_{k, lr,p}}{\partial \beta_{j,o}} =& \, \frac{ \exp_{ mr, p }( \cdot ) \exp_{ mr, o }( \cdot ) \cdot x_j } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ mr, p }( \cdot ) ]^2 } \\ & - \frac{ \exp_{ ml, p }( \cdot ) \exp_{ ml, o }( \cdot ) \cdot x_j } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ ml, p }( \cdot ) ]^2 } \;\;\; \forall \; o \neq p, \; j \in \{ 1, \ldots, K \} \setminus k, k+1 \nonumber \\[1em] \frac{ \partial E_{k, lr, p}}{\partial \delta_{m,p}} =& \, \frac{ \exp_{ ml, p }( \cdot ) ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ ml, o } ( \cdot ) ) \cdot D_{ml} } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ ml, p }( \cdot ) ]^2 } \\ & - \frac{ \exp_{ mr, p }( \cdot ) ( 1 + \sum_{o \in \{ 1, \ldots, P \} \setminus p, p^* } \exp_{ mr, o } ( \cdot ) ) \cdot D_{mr} } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ mr, p }( \cdot ) ]^2 } \nonumber \\[1em] \frac{ \partial E_{k, ml,p}}{\partial \delta_{m,o}} =& \, \frac{ \exp_{mr,p}( \cdot ) \exp_{ mr, o }( \cdot ) \cdot D_{mr} } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ mr, p }( \cdot ) ]^2 } - \frac{ \exp_{ ml, p }( \cdot ) \exp_{ ml, o }( \cdot ) \cdot D_{ml} } {[ 1 + \sum_{p \in \{ 1, \ldots, P \} \setminus p^*} \exp_{ ml, p }( \cdot ) ]^2 } \label{eq:MNLgrad3.4} \end{align} with \begin{align} \exp_{ml} (\cdot) & \equiv \exp \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m D_{ml} \right) \\ \exp_{mr} (\cdot) & \equiv \exp \left( \beta_0 + \sum_{j \in \{ 1, \ldots , K \} \setminus k} \beta_j x_j + \sum_{m=1}^M \delta_m D_{mr} \right) \label{eq:34grad} \end{align} and for the conditional logit as \begin{align} \frac{\partial E_{t,lr}}{\partial \gamma_0} =& \; \frac{\exp_{ml,p}(\cdot) \cdot \sum_{p \in C } \exp_{ml,p} - \exp_{ml,p}(\cdot) \cdot \sum_{p \in C } \exp_{ml,p}} {[\sum_{p \in C } \exp_{ml}(\cdot)]^2} \label{eq:Condgrad3.4a}\\ &- \frac{\exp_{mr,p}(\cdot) \cdot \sum_{p \in C } \exp_{mr,p} - \exp_{mr,p}(\cdot) \cdot \sum_{p \in C } \exp_{mr,p}} {[\sum_{p \in C } \exp_{mr}(\cdot)]^2} \nonumber \\[1em] \frac{\partial E_{t,lr}}{\partial \gamma_j} =& \; \frac{\exp_{ml,p}(\cdot) z_{j,p} \cdot \sum_{p \in C } \exp_{ml,p} - \exp_{ml,p}(\cdot) \cdot \sum_{p \in C } \exp_{ml,p} z_{j,p}} {[\sum_{p \in C } \exp_{ml}(\cdot)]^2} \\ &- \frac{\exp_{mr,p}(\cdot) z_{j,p} \cdot \sum_{p \in C } \exp_{mr,p} - \exp_{mr,p}(\cdot) \cdot \sum_{p \in C } \exp_{mr,p} z_{j,p}} {[\sum_{p \in C } \exp_{mr}(\cdot)]^2} \; \forall \; j \in \{ 1, \ldots , T \} \setminus t \nonumber \\[1em] \frac{\partial E_{t,lr}}{\partial \delta_m} =& \; \frac{\exp_{ml,p}(\cdot) D_{ml} \cdot \sum_{p \in C } \exp_{ml,p} - \exp_{ml,p}(\cdot) \cdot \sum_{p \in C } \exp_{ml,p} D_{ml}} {[\sum_{p \in C } \exp_{ml}(\cdot)]^2} \\ &- \frac{\exp_{mr,p}(\cdot) D_{mr} \cdot \sum_{p \in C } \exp_{mr,p} - \exp_{mr,p}(\cdot) \cdot \sum_{p \in C } \exp_{mr,p} D_{mr}} {[\sum_{p \in C } \exp_{mr}(\cdot)]^2} \; \forall \; m = 1, \ldots, M \label{eq:Condgrad3.4b} \end{align} and for the nested logit as \begin{align} \frac{\partial E_{o,p,lr}}{\partial \beta_0 } =& \; \left( \left( \frac{1}{\lambda_o} - \frac{\sum_{p \in B_o} ( \exp_{o,p,l}(\cdot) \frac{1}{\lambda_o}) } {\sum_{p \in B_o} \exp_{o,p,l}(\cdot)} \right) + \right. \nonumber \\ & \left. \left( \frac{ \sum_{p \in B_o} ( \exp_{ o,p,l} ( \cdot) \frac{1}{\lambda_o} )} { \sum_{p \in B_o} \exp_{o,p,l} ( \cdot)} - \frac{ \sum_{o=1}^O \left( \exp_{o,l} ( \cdot) \frac{ \sum_{p \in B_o} ( \exp{ o,p,l} ( \cdot) 1/\lambda_o )} { \sum_{p \in B_o \exp_{o,p,l} ( \cdot )}} \right)} { \sum_{o=1}^O \exp_{o,l} ( \cdot) } \right) \right) \cdot \pi_{p,l} \pi_{o,l} \nonumber \\ & - \left( \left( \frac{1}{\lambda_o} - \frac{\sum_{p \in B_o} ( \exp_{o,p,r}(\cdot) \frac{1}{\lambda_o}) } {\sum_{p \in B_o} \exp_{o,p,r}(\cdot)} \right) + \right. \nonumber \\ & \left. \left( \frac{ \sum_{p \in B_o} ( \exp_{ o,p,r} ( \cdot) \frac{1}{\lambda_o} )} { \sum_{p \in B_o} \exp_{o,p,r} ( \cdot)} - \frac{ \sum_{o=1}^O \left( \exp_{o,r} ( \cdot) \frac{ \sum_{p \in B_o} ( \exp{ o,p,r} ( \cdot) 1/\lambda_o )} { \sum_{p \in B_o \exp_{o,p,r} ( \cdot )}} \right)} { \sum_{o=1}^O \exp_{o,r} ( \cdot) } \right) \right) \cdot \pi_{p,r} \pi_{o,r} \\[1em] \frac{\partial E_{o,p,lr}}{\partial \beta_j } =& \; \left( \left( \frac{x_{j,p}}{\lambda_o} - \frac{\sum_{p \in B_o} ( \exp_{o,p,l}(\cdot) \frac{ x_{j,p}}{\lambda_o}) } {\sum_{p \in B_o} \exp_{o,p,l}(\cdot)} \right) + \right. \nonumber \\ & \left. \left( \frac{ \sum_{p \in B_o} ( \exp_{ o,p,l} ( \cdot) \frac{x_{p,j}}{\lambda_o} )} { \sum_{p \in B_o} \exp_{o,p,l} ( \cdot)} - \frac{ \sum_{o=1}^O \left( \exp_{o,l} ( \cdot) \frac{ \sum_{p \in B_o} ( \exp{ o,p,l} ( \cdot) x_{j,p}/\lambda_o )} { \sum_{p \in B_o \exp_{o,p,l} ( \cdot )}} \right)} { \sum_{o=1}^O \exp_{o,l} ( \cdot) } \right) \right) \cdot \pi_{p,l} \pi_{o,l} \nonumber \\ & - \left( \left( \frac{x_{j,p}}{\lambda_o} - \frac{\sum_{p \in B_o} ( \exp_{o,p,r}(\cdot) \frac{ x_{j,p}}{\lambda_o}) } {\sum_{p \in B_o} \exp_{o,p,r}(\cdot)} \right) + \right. \nonumber \\ & \left. \left( \frac{ \sum_{p \in B_o} ( \exp_{ o,p,r} ( \cdot) \frac{x_{p,j}}{\lambda_o} )} { \sum_{p \in B_o} \exp_{o,p,r} ( \cdot)} - \frac{ \sum_{o=1}^O \left( \exp_{o,r} ( \cdot) \frac{ \sum_{p \in B_o} ( \exp{ o,p,r} ( \cdot) x_{j,p}/\lambda_o )} { \sum_{p \in B_o \exp_{o,p,r} ( \cdot )}} \right)} { \sum_{o=1}^O \exp_{o,r} ( \cdot) } \right) \right) \cdot \pi_{p,r} \pi_{o,r} \\[1em] \frac{\partial E_{o,p,lr}}{\partial \delta_m } =& \; \left( \left( \frac{D_{m,p,l}}{\lambda_o} - \frac{\sum_{p \in B_o} ( \exp_{o,p,l}(\cdot) \frac{ D_{m,p,l}}{\lambda_o}) } {\sum_{p \in B_o} \exp_{o,p,l}(\cdot)} \right) + \right. \nonumber \\ & \left. \left( \frac{ \sum_{p \in B_o} ( \exp_{ o,p,l} ( \cdot) \frac{D_{m,p,l}}{\lambda_o} )} { \sum_{p \in B_o} \exp_{o,p,l} ( \cdot)} - \frac{ \sum_{o=1}^O \left( \exp_{o,l} ( \cdot) \frac{ \sum_{p \in B_o} ( \exp{ o,p,l} ( \cdot) D_{m,p,l}/\lambda_o )} { \sum_{p \in B_o \exp_{o,p,l} ( \cdot )}} \right)} { \sum_{o=1}^O \exp_{o,l} ( \cdot) } \right) \right) \cdot \pi_{p,l} \pi_{o,l} \nonumber \\ & - \left( \left( \frac{x_{j,p}}{\lambda_o} - \frac{\sum_{p \in B_o} ( \exp_{o,p,r}(\cdot) \frac{ D_{m,p,r}}{\lambda_o}) } {\sum_{p \in B_o} \exp_{o,p,r}(\cdot)} \right) + \right. \nonumber \\ & \left. \left( \frac{ \sum_{p \in B_o} ( \exp_{ o,p,r} ( \cdot) \frac{D_{m,p,r}}{\lambda_o} )} { \sum_{p \in B_o} \exp_{o,p,r} ( \cdot)} - \frac{ \sum_{o=1}^O \left( \exp_{o,r} ( \cdot) \frac{ \sum_{p \in B_o} ( \exp{ o,p,r} ( \cdot) D_{m,p,r}/\lambda_o )} { \sum_{p \in B_o \exp_{o,p,r} ( \cdot )}} \right)} { \sum_{o=1}^O \exp_{o,r} ( \cdot) } \right) \right) \cdot \pi_{p,r} \pi_{o,r} \end{align} with \begin{align} \exp_{p,o,l}(\cdot) &\equiv \exp \left( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \sum_{m=1}^M \delta_m/\lambda_o D_{m,p,l} \right) \\[1em] \exp_{p,o,r}(\cdot) &\equiv \exp \left( \beta_0/\lambda_o + \sum_{j \in \{1, \ldots, K \} \setminus k } \beta_j/\lambda_o x_{j,p} + \sum_{m=1}^M \delta_m/\lambda_o D_{m,p,l} \right) \\[1em] \exp_{l/r,o} (\cdot) &\equiv \exp \left( \gamma_0 + \sum_{t=1}^T \gamma_t z_{t,o} + \lambda_o IV_o \right) \end{align} The following function calculates the effect $E_{k,lr}$ and its standard error in accordance with equations~(\ref{eq:logitEffGr}), and (\ref{eq:34Delta}) to (\ref{eq:34grad}) for the binary logit, effect $E_{k,lr,p}$ and its standard error in accordance with equations~(\ref{eq:MNlogitEffGr}) and (\ref{eq:MNLGrad1}) to (\ref{eq:MNLgrad3.4}) for the multi-nomial logit, and effect $E_{t,lr,p}$ and its standard error in accordance with equations~(\ref{eq:CondlogitEffGr}) and (\ref{eq:Condgrad3.4a}) to (\ref{eq:Condgrad3.4b}): << >>= logEffGr <- function( allCoef, allXVal, xPos, Group, yCat = NA, allCoefSE = rep( NA, length( allCoef ) ), method = "binary" ){ if( method == "binary" ){ nCoef <- length( allCoef ) xCoef <- allCoef[ xPos ] xShares <- allXVal[ xPos ] } else if( method == "MNL" ){ nCoef <- length( allCoef ) mCoef <- matrix( allCoef, nrow = length( allXVal ) ) NCoef <- dim( mCoef )[2] pCoef <- dim( mCoef )[1] xCoef <- mCoef[ xPos, ] xShares <- allXVal[ xPos ] } else{ nCoef <- length( allCoef ) xCoef <- allCoef[ xPos ] mXVal <- matrix( allXVal, nrow = nCoef ) pCoef <- dim( mXVal )[2] xShares <- mXVal[ xPos, ] } if( method == "binary" || method == "MNL" ){ if( sum( xShares ) > 1 ){ stop( "the shares in argument 'xShares' sum up to a value larger than 1" ) } } else{ for( p in 1:pCoef ){ if( sum( xShares[ , p ] ) > 1 ){ stop( "the shares in argument 'xShares' sum up to a value larger than 1" ) } } } if( method == "binary" ){ if( length( xCoef ) != length( xShares ) ){ stop( "arguments 'xCoef' and 'xShares' must have the same length" ) } if( length( xCoef ) != length( Group ) ){ stop( "arguments 'xCoef' and 'Group' must have the same length" ) } } else if( method == "MNL" ){ if( dim( xCoef )[1] != length( xShares ) ){ stop( "arguments 'xCoef' and 'xShares' must have the same length" ) } if( dim( xCoef )[1] != length( Group ) ){ stop( "arguments 'xCoef' and 'Group' must have the same length" ) } } else{ if( length( xCoef ) != dim( xShares )[1] ){ stop( "arguments 'xCoef' and 'xShares' must have the same length" ) } if( length( xCoef ) != length( Group ) ){ stop( "arguments 'xCoef' and 'Group' must have the same length" ) } } if( !all( Group %in% c( -1, 0, 1 ) ) ){ stop( "all elements of argument 'Group' must be -1, 0, or 1" ) } if( method == "binary" ){ # D_mr DRef <- sum( xCoef[ Group == -1 ] * xShares[ Group == -1 ]) / sum( xShares[ Group == -1 ] ) XBetaRef <- sum( allCoef[ -xPos ] * allXVal[ -xPos ]) + DRef # D_ml DEffect <- sum( xCoef[ Group == 1 ] * xShares[ Group == 1 ]) / sum( xShares[ Group == 1 ] ) XBetaEffect <- sum( allCoef[ -xPos ] * allXVal[ -xPos ]) + DEffect # effect effeG <- exp( XBetaEffect )/( 1 + exp( XBetaEffect ) ) - exp( XBetaRef )/( 1 + exp( XBetaEffect ) ) } else if( method == "MNL" ){ # D_mr DRef <- colSums( xCoef[ Group == -1, , drop = FALSE ] * xShares[ Group == -1 ] )/ sum( xShares[ Group == -1 ] ) XBetaRef <- colSums( mCoef[ -xPos, , drop = FALSE ] * allXVal[ -xPos ]) + DRef # D_ml DEffect <- colSums( xCoef[ Group == 1, , drop = FALSE ] * xShares[ Group == 1 ] )/ sum( xShares[ Group == 1 ] ) XBetaEffect <- colSums( mCoef[ -xPos, , drop = FALSE ] * allXVal[ -xPos ]) + DEffect # effect effeG <- exp( XBetaEffect[ yCat ] )/( 1 + sum( exp( XBetaEffect ) ) ) - exp( XBetaRef[ yCat ] )/( 1 + sum( exp( XBetaRef ) ) ) } else{ # D_mr DRef <- colSums( xCoef[ Group == -1 ] * xShares[ Group == -1, , drop = FALSE ] )/ sum( xShares[ Group == -1, , drop = FALSE ] ) XBetaRef <- colSums( allCoef[ -xPos ] * mXVal[ -xPos, , drop = FALSE ] ) + DRef # D_ml DEffect <- colSums( xCoef[ Group == 1 ] * xShares[ Group == 1, , drop = FALSE ] )/ sum( xShares[ Group == 1, , drop = FALSE ] ) XBetaEffect <- colSums( allCoef[ -xPos ] * mXVal[ -xPos, , drop = FALSE ] ) + DEffect # effect effeG <- exp( XBetaEffect[ yCat ] )/( sum( exp( XBetaEffect ) ) ) - exp( XBetaRef[ yCat ] )/( sum( exp( XBetaRef ) ) ) } # partial derivative of E_{k,ml} w.r.t. all estimated coefficients if( method == "binary" ){ derivCoef <- rep( NA, nCoef ) derivCoef[ -xPos ] = ( exp( XBetaEffect )/( 1 + exp( XBetaEffect ))^2 - exp( XBetaRef )/( 1 + exp( XBetaRef ))^2 ) * allXVal[ -xPos ] derivCoef[ xPos ] = exp( XBetaEffect )/( 1 + exp( XBetaEffect))^2 * DEffect - exp( XBetaRef )/( 1 + exp( XBetaRef ))^2 * DRef } else if( method == "MNL" ){ derivCoef <- matrix( NA, nrow=pCoef, ncol=NCoef ) for( p in 1:NCoef ){ if( p == yCat ){ derivCoef[ -xPos, p ] <- ( exp( XBetaEffect[ p ] ) * ( 1 + sum( exp( XBetaEffect[ -yCat ] ) ) )/ ( 1 + sum( exp( XBetaEffect ) ) )^2 - exp( XBetaRef[ p ] ) * ( 1 + sum( exp( XBetaRef[ -yCat ] ) ) )/ ( 1 + sum( exp( XBetaRef ) ) )^2 ) * allXVal[ -xPos ] derivCoef[ xPos, p ] <- ( exp( XBetaEffect[ p ] ) * ( 1 + sum( exp( XBetaEffect[ -yCat ] ) ) )/ ( 1 + sum( exp( XBetaEffect ) ) )^2 ) * DEffect - ( exp( XBetaRef[ p ] ) * ( 1 + sum( exp( XBetaRef[ -yCat ] ) ) )/ ( 1 + sum( exp( XBetaRef ) ) )^2 ) * DRef } else{ derivCoef[ -xPos, p ] <- ( ( exp( XBetaRef[ yCat ] ) * exp( XBetaRef[ p ] ) )/ ( 1 + sum( exp( XBetaRef ) ) )^2 - ( exp( XBetaEffect[ yCat ] ) * exp( XBetaEffect[ p ] ) )/ ( 1 + sum( exp( XBetaEffect ) ) )^2 ) * allXVal[ -xPos ] derivCoef[ xPos, p ] <- ( ( exp( XBetaRef[ yCat ] ) * exp( XBetaRef[ p ] ) )/ ( 1 + sum( exp( XBetaRef ) ) )^2 ) * DRef - ( ( exp( XBetaEffect[ yCat ] ) * exp( XBetaEffect[ p ] ) )/ ( 1 + sum( exp( XBetaEffect ) ) )^2 ) * DEffect } } derivCoef <- c( derivCoef ) } else{ derivCoef <- rep( NA, nCoef ) derivCoef[ -xPos ] = ( exp( XBetaEffect[ yCat ] ) * mXVal[ -xPos, yCat ] * sum( exp( XBetaEffect ) ) - exp( XBetaEffect[ yCat ] ) * sum( exp( XBetaEffect ) * mXVal[ -xPos, ] ) )/ ( sum( exp( XBetaEffect ) ) )^2 - ( exp( XBetaRef[ yCat ] ) * mXVal[ -xPos, yCat ] * sum( exp( XBetaRef ) ) - exp( XBetaRef[ yCat ] ) * sum( exp( XBetaRef ) * mXVal[ -xPos, ] ) )/ ( sum( exp( XBetaRef ) ) )^2 derivCoef[ xPos ] = ( exp( XBetaEffect[ yCat ] ) * DEffect[ yCat ] * sum( exp( XBetaEffect ) ) - exp( XBetaEffect[ yCat ] ) * sum( exp( XBetaEffect ) * DEffect[ yCat ] ) )/ ( sum( exp( XBetaEffect ) ) )^2 - ( exp( XBetaRef[ yCat ] ) * DRef[ yCat ] * sum( exp( XBetaRef ) ) - exp( XBetaRef[ yCat ] ) * sum( exp( XBetaRef ) * DRef[ yCat ] ) )/ ( sum( exp( XBetaRef ) ) )^2 } # variance covariance of the coefficients (covariances set to zero) vcovCoef <- diag( allCoefSE^2 ) # approximate standard error of the effect effeGSE <- drop( sqrt( t( derivCoef ) %*% vcovCoef %*% derivCoef ) ) result <- c( effect = effeG, stdEr = effeGSE ) return( result ) } @ where argument \code{allCoef} $= (\beta_0, \ldots, \beta_{k-1}, \beta_{k+1}, \ldots, \beta_K, \delta_1, \ldots, \delta_M )^\top$ are the coefficients of the binary or conditional logit regression or \code{allCoef} $= ( \beta_{0,1}, \ldots, \beta_{k-1,1}, \beta_{k+1,1}, \ldots, \beta_{K,1}, \delta_{1,1}, \ldots, \delta_{M,1}, \beta_{0,P}, \ldots, \beta_{k-1,P}, \beta_{k+1,P}, \ldots, \beta_{K,P}, \delta_{1,P}, \ldots, \delta_{M,P} )^\top$ a vector of all the $P$ sets of coefficients from the multi-nomial logit regression which does \textbf{not} include any values for the reference category, \textbf{with the coefficient of the reference group of the $k$th variable set to zero}; argument \code{allXVal} $= (x_1, \ldots, x_{k-1}, x_{k+1}, \ldots, x_K, s_1, \ldots, s_M )^\top$ are the mean values of the respective explanatory variables and the shares of each of the categories of the $k$th variable \textbf{which includes the share of the reference group of the $k$th variable} or \code{allXVal} $= ( z_{0,1}, \ldots, z_{t-1,1}, z_{t+1,1}, \ldots, z_{T,1}, s_{1,1}, \ldots, s_{M,1}, \ldots, z_{0,P}, \ldots, z_{t-1,P}, z_{t+1,P}, \ldots, z_{T,P}, s_{1,P}, \ldots, s_{M,P} )^\top$ are the mean values of the $P$ sets of explanatory variables and the $P$ sets of shares of the $t$th variable \textbf{which includes the share of the reference group of the $t$th variable}; argument \code{xPos} is a vector of at least two elements indicating the position of the $k$th or $t$th variable in arguments \code{allCoef} and \code{allXVal}; argument \code{Group} is a vector of at least two elements consisting of $-1$s, $0$s, and $1$s, where a $-1$ indicates that the category belongs to the new reference category, a $1$ indicates that the category belongs to the new category of interest, and a zero indicates that the category belongs to neither; argument \code{allCoefSE} $=( \se( \beta_0), \ldots, \se( \beta_{k-1}), \se( \beta_{k+1}), \ldots, \se( \beta_K), \se(\delta_1), \ldots, \se(\delta_M ))^\top$ are the standard errors of all coefficients of the binary or conditional logit regression or \code{allCoefSE} $=( \se(\beta_{0,1}), \ldots, \se(\beta_{k-1,1}), \se(\beta_{k+1,1}), \ldots, \se(\beta_{K,1}), \se(\delta_{1,1}), \ldots, \se(\delta_{M,1}), \se(\beta_{0,P}), \ldots, \se(\beta_{k-1,P}), \se(\beta_{k+1,P}), \ldots, \se(\beta_{K,P}), \se(\delta_{1,P}), \ldots, \se(\delta_{M,P} ) )^\top$ of the multi-nomial logit, respectively, \textbf{where the standard error for the reference group is included as zero}; argument \code{method = "binary", "MNL", "CondL" } identifies the estimator chosen, where \code{"binary"} indicates the binary logit estimator, \code{"MNL"} indicates the multi-nomial logit estimator, and \code{"CondL"} indicates the conditional logit; \code{yCat} is only relevant in connection with method \code{"MNL"} or \code{"CondL"} and indicates the $p$th output category of the multi-nomial logit that is of interest. Elements of arguments \code{allCoef}, \code{allXVal}, \code{Group}, and \code{allCoefSE} that belong to a category that is neither in the reference category nor in the category of interest, i.e., categories~$m$ with $D_{mr} = D_{ml} = 0$, can be omitted but the ommission must be consistent for all four arguments. << >>= # Example eff10a <- logEffGr( allCoef = c( 0.28, 0.003, 0.0075, 0, -0.034, -0.005, 0.89, -1.2 ), allXVal = c( 1, 0.1, 0.3, 0.25, 0.15, 0.2, 2.34, 10.8 ), xPos = c( 2:6 ), Group = c( 0, -1, -1, 1, 1 ) ) eff10a eff10b <- logEffGr( allCoef = c( 0.28, 0.003, 0.0075, 0, -0.034, -0.005, 0.89, -1.2 ), allXVal = c( 1, 0.1, 0.3, 0.25, 0.15, 0.2, 2.34, 10.8 ), xPos = c( 2:6 ), Group = c( 0, -1, -1, 1, 1 ), allCoefSE = c( 0.03, 0.0001, 0.005, 0, 0.01, 0.004, 0.05, 0.8 ) ) eff10b # Example eff11a <- logEffGr( allCoef = c( 0.28, 0.003, 0.0075, 0, -0.034, -0.005, 0.89, 0, 0.005, 0.06, 1.7, 0 ), allXVal = c( 1, 0.5, 0.3, 0.2 ), xPos = c( 2:4 ), Group = c( -1, -1, 1 ), yCat = 2, method = "MNL" ) eff11a eff11b <- logEffGr( allCoef = c( 0.28, 0.003, 0.0075, 0, -0.034, -0.005, 0.89, 0, 0.005, 0.06, 1.7, 0 ), allXVal = c( 1, 0.5, 0.3, 0.2 ), xPos = c( 2:4 ), Group = c( -1, -1, 1 ), yCat = 2, method = "MNL", allCoefSE = c( 0.03, 0.0001, 0.005, 0, 0.01, 0.004, 0.05, 0, 0.004, 0.5, 0.0078, 0 ) ) eff11b # Example eff12a <- logEffGr( allCoef = c( 0.28, 0.003, 0.0075, 0 ), allXVal = c( 1, 0.5, 0.3, 0.2, 1, 0.4, 0.4, 0.1 ), xPos = c( 2:4 ), Group = c( -1, -1, 1 ), yCat = 2, method = "CondL" ) eff12a eff12b <- logEffGr( allCoef = c( 0.28, 0.003, 0.0075, 0 ), allXVal = c( 1, 0.5, 0.3, 0.2, 1, 0.4, 0.4, 0.1 ), xPos = c( 2:4 ), allCoefSE = c( 0.03, 0.0001, 0.005, 0 ), Group = c( -1, -1, 1 ), yCat = 2, method = "CondL" ) eff12b @ \clearpage \appendix \noindent{\Huge\textbf{APPENDIX}} \section{Additional helper functions} \label{sec:helperFunctions} The following code defines a helper function that checks argument \code{xPos} of some functions that are defined above: <>= @ The following code defines a helper function that checks whether $\beta_0 + \sum_{j=1}^K \beta_j x_j$ has a plausible value: <>= @ \end{document}