## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ## ----------------------------------------------------------------------------- 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 ) } ## ----elaIntBounds------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ## ----------------------------------------------------------------------------- EXSquared <- function( lowerBound, upperBound ) { result <- ( upperBound^3 - lowerBound^3 )/( 3 * ( upperBound - lowerBound ) ) return( result ) } ## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ## ----checkXPos,echo=FALSE----------------------------------------------------- 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 ) } } } ## ----checkXBeta,echo=FALSE---------------------------------------------------- checkXBeta <- function( xBeta ) { if( any( abs( xBeta ) > 3.5 ) ) { warning( "At least one x'beta has an implausible value: ", paste( xBeta, collapse = ", " ) ) } } ## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ) ) ) ## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ## ----------------------------------------------------------------------------- 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 ) } ## ----------------------------------------------------------------------------- # 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 ## ----checkXPos,eval=FALSE----------------------------------------------------- # 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 ) # } # } # } ## ----checkXBeta,eval=FALSE---------------------------------------------------- # checkXBeta <- function( xBeta ) { # if( any( abs( xBeta ) > 3.5 ) ) { # warning( "At least one x'beta has an implausible value: ", # paste( xBeta, collapse = ", " ) ) # } # }