# R/urbinEla.R In urbin: Unifying Estimation Results with Binary Dependent Variables

```urbinEla <- function( allCoef, allXVal, xPos, model,
allCoefVcov = NULL, seSimplify = !is.matrix( allCoefVcov ),
xMeanSd = NULL, iPos = 1, yCat = NULL ){

# check argument seSimplify
if( length( seSimplify ) != 1 || !is.logical( seSimplify ) ) {
stop( "argument 'seSimplify' must be TRUE or FALSE" )
}
if( model == "lpm" && !is.null( allCoefVcov ) &&
( seSimplify == is.matrix( allCoefVcov ) ) ) {
warning( "argument 'seSimplify' is ignored in linear probability models" )
} else if( seSimplify && is.matrix( allCoefVcov ) ) {
warning( "if the (full) covariance matrix is provided",
" via argument 'allCoefCov'",
" it is NOT recommended to set argument 'seSimplify' to TRUE" )
} else if( !seSimplify && !is.null( allCoefVcov) &&
!is.matrix( allCoefVcov ) ) {
warning( "the returned standard error is likely very imprecise;",
" you can provide the full covariance matrix",
" via argument 'allCoefVcov'",
" or do NOT set argument 'seSimplify' to FALSE",
" to obtain a more precise standard error" )
}
# number of coefficients
nCoef <- length( allCoef )
# number of explanatory variables
nXVal <- length( allXVal )
# check allXVal and allCoef
if( model %in% c( "lpm", "probit", "oprobit", "logit" ) ){
# LPM model: allXVal can be a scalar even if there is a quadratic term
if( model == "lpm" && length( xPos ) == 2 && length( allXVal ) == 1 ){
allXVal <- c( allXVal, allXVal^2 )
nXVal <- length( allXVal )
}
if( nXVal != nCoef ) {
stop( "arguments 'allCoef' and 'allXVal' must have the same length" )
}
} else if( model == "mlogit" ){
# number of alternative categories of the dependent variable
nYCat <- round( nCoef / nXVal )
if( nCoef != nXVal * nYCat ) {
stop( "length of argument 'allCoef' must be a multiple",
" of the length of argument 'allXVal'" )
}
# create matrix of coefficients
mCoef <- matrix( allCoef, nrow = nXVal, ncol = nYCat )
# add column for coefficients of the reference category
mCoef <- cbind( mCoef, 0 )
} else {
stop( "argument 'model' specifies an unknown type of model" )
}
# check argument yCat
if( model == "mlogit" ) {
checkYCat( yCat, nYCat, maxLength = nYCat + 1 )
yCat[ yCat == 0 ] <- nYCat + 1
} else if( !is.null( yCat ) ) {
warning( "argument 'yCat' is ignored" )
}
# Check position vector
checkXPos( xPos, minLength = 1, maxLength = 2, minVal = 1,
maxVal = ifelse( model == "mlogit", nXVal, nCoef ) )
# check position of the intercept
checkIPos( iPos, xPos, allXVal, model )
# check and prepare allCoefVcov
allCoefVcov <- prepareVcov( allCoefVcov, length( allCoef ), xPos, xMeanSd,
nXVal = nXVal, iPos = iPos, pCall = match.call() )
# Identify coefficients of interest (kth/tth covariate)
if( length( xPos ) == 2 ){
if( model %in% c( "lpm", "probit", "oprobit", "logit" ) ){
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( model == "mlogit" ){
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 {
stop( "argument 'model' specifies an unknown type of model" )
}
} else if( length( xPos ) == 1 ) {
if( model %in% c( "lpm", "probit", "oprobit", "logit" ) ){
xCoef <- c( allCoef[ xPos ], 0 )
} else if( model == "mlogit" ){
xCoef <- rbind( mCoef[ xPos, ], 0 )
} else {
stop( "argument 'model' specifies an unknown type of model" )
}
} else {
stop( "argument 'xPos' must be a scalar or a vector with two elements" )
}
# prepare calculation of semi-elasticity
if( model == "lpm" ) {
xVal <- allXVal[ xPos[ 1 ] ]
semEla <- ( xCoef[1] + 2 * xCoef[2] * xVal ) * xVal
} else if( model %in% c( "probit", "oprobit" ) ){
xVal <- allXVal[ xPos[ 1 ] ]
xBeta <- sum( allCoef * allXVal )
checkXBeta( xBeta )
dfun <- dnorm( xBeta )
semEla <- ( xCoef[ 1 ] + 2 * xCoef[ 2 ] * xVal ) * xVal * dfun
} else if( model == "logit" ){
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( model == "mlogit" ){
xVal <- allXVal[ xPos[1] ]
xCoefLinQuad <- xCoef[ 1, ] + 2 * xCoef[ 2, ] * xVal
xBeta <- allXVal %*% mCoef
checkXBeta( xBeta )
pfun <- exp( xBeta ) / sum( exp( xBeta ) )
semElaCat <- xVal * pfun *
semEla <- sum( semElaCat[ yCat ] )
} else {
stop( "argument 'model' specifies an unknown type of model" )
}
# partial derivatives of semi-elasticities wrt coefficients
if( model == "lpm" ){
derivCoef <- rep( 0, nCoef )
derivCoef[ xPos[1] ] <- xVal
if( length( xPos ) == 2 ) {
derivCoef[ xPos[2] ] <- 2 * xVal^2
}
} else if( model %in% c( "probit", "oprobit" ) ){
if( seSimplify ) {
derivCoef <- rep( 0, length( allCoef ) )
} else {
derivCoef <- - semEla * xBeta * allXVal
}
derivCoef[ xPos[1] ] <- derivCoef[ xPos[1] ] + dnorm( xBeta ) * xVal
if( length( xPos ) == 2 ) {
derivCoef[ xPos[2] ] <- derivCoef[ xPos[2] ] + dnorm( xBeta ) * 2 * xVal^2
}
} else if( model == "logit" ){
if( seSimplify ) {
derivCoef <- rep( 0, length( allCoef ) )
} else {
derivCoef <- allXVal * semEla *
( 1 - 2 * exp( xBeta ) / ( 1 + exp( xBeta ) ) )
}
derivCoef[ xPos[1] ] <- derivCoef[ xPos[1] ] + dfun * xVal
if( length( xPos ) == 2 ) {
derivCoef[ xPos[2] ] <- derivCoef[ xPos[2] ] + dfun * 2 * xVal^2
}
} else if( model == "mlogit" ){
derivCoef <- rep( 0, length( allCoef ) )
for( p in 1:nYCat ) {
coefNoYCat <- ( 1 + (p-1)*nXVal ):( p * nXVal )
for( yCati in yCat ) {
derivCoef[ coefNoYCat ][ xPos[1] ] <-
derivCoef[ coefNoYCat ][ xPos[1] ] -
pfun[ yCati ] * xVal * pfun[ p ] +
( p == yCati ) * pfun[ yCati ] * xVal
if( length( xPos ) == 2 ) {
derivCoef[ coefNoYCat ][ xPos[2] ] <-
derivCoef[ coefNoYCat ][ xPos[2] ] -
pfun[ yCati ] * 2 * xVal^2 * pfun[ p ] +
( p == yCati ) * pfun[ yCati ] * 2 * xVal^2
}
if( !seSimplify ) {
derivCoef[ coefNoYCat ][ -xPos ] <-
derivCoef[ coefNoYCat ][ -xPos ] +
( - pfun[ yCati ] * semElaCat[ p ] -
pfun[ p ] * semElaCat[ yCati ] +
( p == yCati ) * semElaCat[ p ] ) *
allXVal[ -xPos ]
derivCoef[ coefNoYCat ][ xPos[1] ] <-
derivCoef[ coefNoYCat ][ xPos[1] ] +
( - pfun[ yCati ] * semElaCat[ p ] -
pfun[ p ] * semElaCat[ yCati ] +
( p == yCati ) * semElaCat[ p ] ) *
xVal
if( length( xPos ) == 2 ) {
derivCoef[ coefNoYCat ][ xPos[2] ] <-
derivCoef[ coefNoYCat ][ xPos[2] ] +
( - pfun[ yCati ] * semElaCat[ p ] -
pfun[ p ] * semElaCat[ yCati ] +
( p == yCati ) * semElaCat[ p ] ) *
xVal^2
}
}
}
}
} else {
stop( "argument 'model' specifies an unknown type of model" )
}

# approximate standard error of the semi-elasticity
semElaSE <- drop( sqrt( t( derivCoef ) %*% allCoefVcov %*% derivCoef ) )

# create object that will be returned
result <- list()
result\$call <- match.call()
result\$allCoefVcov <- allCoefVcov
result\$derivCoef <- derivCoef
result\$semEla <- unname( semEla )
result\$stdEr <- semElaSE
class( result ) <- "urbin"
return( result )
}
```

## Try the urbin package in your browser

Any scripts or data that you put into this service are public.

urbin documentation built on May 2, 2019, 5:23 p.m.