R Under development (unstable) (2024-09-30 r87204 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library( "urbin" ) > maxLikLoaded <- require( "maxLik" ) Loading required package: maxLik Loading required package: miscTools Please cite the 'maxLik' package as: Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1. If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site: https://r-forge.r-project.org/projects/maxlik/ > mfxLoaded <- require( "mfx" ) Loading required package: mfx Loading required package: sandwich Loading required package: lmtest Loading required package: zoo Attaching package: 'zoo' The following objects are masked from 'package:base': as.Date, as.Date.numeric Loading required package: MASS Loading required package: betareg > if( !require( "sampleSelection" ) ) { + q( save = "no" ) + } Loading required package: sampleSelection > options( digits = 4 ) > > # load data set > data( "Mroz87", package = "sampleSelection" ) > > # create dummy variable for kids > Mroz87$kids <- as.numeric( Mroz87$kids5 > 0 | Mroz87$kids618 > 0 ) > > ### linear in age > estProbitLin <- glm( lfp ~ kids + age + educ, + family = binomial(link = "probit"), + data = Mroz87 ) > summary( estProbitLin ) Call: glm(formula = lfp ~ kids + age + educ, family = binomial(link = "probit"), data = Mroz87) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.07161 0.45249 -0.16 0.8743 kids -0.31563 0.12223 -2.58 0.0098 ** age -0.01916 0.00698 -2.74 0.0061 ** educ 0.10464 0.02124 4.93 8.3e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1029.75 on 752 degrees of freedom Residual deviance: 993.43 on 749 degrees of freedom AIC: 1001 Number of Fisher Scoring iterations: 4 > # mean values of the explanatory variables > xMeanLin <- c( 1, colMeans( Mroz87[ , c( "kids", "age", "educ" ) ] ) ) > # semi-elasticity of age without standard errors > urbinEla( coef( estProbitLin ), xMeanLin, xPos = 3, model = "probit" ) semEla stdEr -0.3199 NA > # semi-elasticity of age based on numerical derivation > 100 * ( predict( estProbitLin, + newdata = as.data.frame( t( xMeanLin * c( 1, 1, 1.005, 1 ) ) ), + type = "response" ) - + predict( estProbitLin, + newdata = as.data.frame( t( xMeanLin * c( 1, 1, 0.995, 1 ) ) ), + type = "response" ) ) 1 -0.3199 > # partial derivatives of the semi-elasticity wrt the coefficients > urbinEla( coef( estProbitLin ), xMeanLin, 3, model = "probit", + seSimplify = FALSE )$derivCoef kids age educ 0.05743 0.03996 19.14172 0.70559 > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients > if( maxLikLoaded ) { + print( numericGradient( function( x, ... ){ urbinEla( x, ... )$semEla }, + t0 = coef( estProbitLin ), + allXVal = xMeanLin, xPos = 3, model = "probit" ) ) + } (Intercept) kids age educ [1,] 0.05743 0.03996 19.14 0.7056 > # simplified partial derivatives of the semi-elasticity wrt the coefficients > urbinEla( coef( estProbitLin ), xMeanLin, 3, model = "probit", + seSimplify = TRUE )$derivCoef [1] 0.0 0.0 16.7 0.0 > # semi-elasticity of age with standard errors (full covariance matrix) > urbinEla( coef( estProbitLin ), xMeanLin, 3, model = "probit", + vcov( estProbitLin ) ) semEla stdEr -0.3199 0.1166 > # semi-elasticity of age with standard errors (only standard errors) > urbinEla( coef( estProbitLin ), xMeanLin, 3, model = "probit", + sqrt( diag( vcov( estProbitLin ) ) ), seSimplify = FALSE ) semEla stdEr -0.3199 0.1371 Warning message: In urbinEla(coef(estProbitLin), xMeanLin, 3, model = "probit", sqrt(diag(vcov(estProbitLin))), : 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 > # semi-elasticity of age with standard errors (only standard errors, simplified) > urbinEla( coef( estProbitLin ), xMeanLin, 3, model = "probit", + sqrt( diag( vcov( estProbitLin ) ) ) ) semEla stdEr -0.3199 0.1166 > # semi-elasticity of age based on partial derivative calculated by the mfx package > if( mfxLoaded ) { + estProbitLinMfx <- probitmfx( lfp ~ kids + age + educ, data = Mroz87 ) + print( estProbitLinMfx$mfxest[ "age", 1:2 ] * xMeanLin[ "age" ] ) + } dF/dx Std. Err. -0.3199 0.1166 > if( mfxLoaded ) { + print( urbinEla( estProbitLinMfx$mfxest[ "age", 1 ], xMeanLin["age"], 1, iPos = 0, + model = "lpm", estProbitLinMfx$mfxest[ "age", 2 ] ) ) + } semEla stdEr -0.3199 0.1166 > if( mfxLoaded ) { + print( urbinEla( estProbitLinMfx$mfxest[ , 1 ], xMeanLin[-1], 2, iPos = 0, + model = "lpm", estProbitLinMfx$mfxest[ , 2 ] ) ) + } semEla stdEr -0.3199 0.1166 > > > ### quadratic in age > estProbitQuad <- glm( lfp ~ kids + age + I(age^2) + educ, + family = binomial(link = "probit"), + data = Mroz87 ) > summary( estProbitQuad ) Call: glm(formula = lfp ~ kids + age + I(age^2) + educ, family = binomial(link = "probit"), data = Mroz87) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.328938 1.394393 -3.10 0.00191 ** kids -0.438760 0.129528 -3.39 0.00071 *** age 0.192165 0.065880 2.92 0.00354 ** I(age^2) -0.002493 0.000773 -3.22 0.00126 ** educ 0.107041 0.021411 5.00 5.8e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1029.75 on 752 degrees of freedom Residual deviance: 982.88 on 748 degrees of freedom AIC: 992.9 Number of Fisher Scoring iterations: 4 > # mean values of the explanatory variables > xMeanQuad <- c( xMeanLin[ 1:3], xMeanLin[3]^2, xMeanLin[4] ) > # semi-elasticity of age without standard errors > urbinEla( coef( estProbitQuad ), xMeanQuad, c( 3, 4 ), model = "probit" ) semEla stdEr -0.3193 NA > # semi-elasticity of age based on numerical derivation > 100 * ( predict( estProbitQuad, + newdata = as.data.frame( t( xMeanQuad * c( 1, 1, 1.005, 1.005^2, 1 ) ) ), + type = "response" ) - + predict( estProbitQuad, + newdata = as.data.frame( t( xMeanQuad * c( 1, 1, 0.995, 0.995^2, 1 ) ) ), + type = "response" ) ) 1 -0.3193 > # partial derivatives of the semi-elasticity wrt the coefficients > urbinEla( coef( estProbitQuad ), xMeanQuad, c( 3, 4 ), model = "probit", + seSimplify = FALSE )$derivCoef kids age age educ 1.097e-01 7.634e-02 2.066e+01 1.560e+03 1.348e+00 > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients > if( maxLikLoaded ) { + print( numericGradient( function( x, ... ){ urbinEla( x, ... )$semEla }, + t0 = coef( estProbitQuad ), + allXVal = xMeanQuad, xPos = c( 3, 4 ), model = "probit" ) ) + } (Intercept) kids age I(age^2) educ [1,] 0.1097 0.07634 20.66 1560 1.348 > # simplified partial derivatives of the semi-elasticity wrt the coefficients > urbinEla( coef( estProbitQuad ), xMeanQuad, c( 3, 4 ), model = "probit", + seSimplify = TRUE )$derivCoef [1] 0 0 16 1361 0 > # semi-elasticity of age with standard errors (full covariance matrix) > urbinEla( coef( estProbitQuad ), xMeanQuad, c( 3, 4 ), model = "probit", + vcov( estProbitQuad ) ) semEla stdEr -0.3193 0.1120 > # semi-elasticity of age with standard errors (only standard errors) > urbinEla( coef( estProbitQuad ), xMeanQuad, c( 3, 4 ), model = "probit", + sqrt( diag( vcov( estProbitQuad ) ) ), seSimplify = FALSE ) semEla stdEr -0.3193 1.8254 Warning messages: 1: In urbinEla(coef(estProbitQuad), xMeanQuad, c(3, 4), model = "probit", : 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 2: In urbinEla(allCoef = coef(estProbitQuad), allXVal = xMeanQuad, xPos = c(3, 4), model = "probit", allCoefVcov = sqrt(diag(vcov(estProbitQuad))), seSimplify = FALSE) : the returned standard error is likely largely upward biased and, thus, in most cases meaningless; you can provide the full covariance matrix via argument 'allCoefVcov' to avoid this bias or use argument 'xMeanSd' to substantially reduce this bias > # semi-elasticity of age with standard errors (only standard errors, simplified) > urbinEla( coef( estProbitQuad ), xMeanQuad, c( 3, 4 ), model = "probit", + sqrt( diag( vcov( estProbitQuad ) ) ) ) semEla stdEr -0.3193 1.4895 Warning message: In urbinEla(allCoef = coef(estProbitQuad), allXVal = xMeanQuad, xPos = c(3, 4), model = "probit", allCoefVcov = sqrt(diag(vcov(estProbitQuad)))) : the returned standard error is likely largely upward biased and, thus, in most cases meaningless; you can provide the full covariance matrix via argument 'allCoefVcov' to avoid this bias or use argument 'xMeanSd' to substantially reduce this bias > # semi-elasticity of age with standard errors (only standard errors, xMeanSd) > urbinEla( coef( estProbitQuad ), xMeanQuad, c( 3, 4 ), model = "probit", + sqrt( diag( vcov( estProbitQuad ) ) ), + xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ), + seSimplify = FALSE ) semEla stdEr -0.3193 0.1105 Warning message: In urbinEla(coef(estProbitQuad), xMeanQuad, c(3, 4), model = "probit", : 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 > # semi-elasticity of age with standard errors (only standard errors, xMeanSd, simplified) > urbinEla( coef( estProbitQuad ), xMeanQuad, c( 3, 4 ), model = "probit", + sqrt( diag( vcov( estProbitQuad ) ) ), + xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) ) semEla stdEr -0.3193 0.1352 > # semi-elasticity of age based on partial derivatives calculated by the mfx package > # (differs from the above, because mean(age)^2 is not the same as mean(age^2)) > if( mfxLoaded ) { + estProbitQuadMfx <- probitmfx( lfp ~ kids + age + I(age^2) + educ, data = Mroz87 ) + print( estProbitQuadMfx$mfxest[ "age", 1:2 ] * xMeanQuad[ "age" ] + + 2 * estProbitQuadMfx$mfxest[ "I(age^2)", 1:2 ] * xMeanQuad[ "age" ]^2 ) + } dF/dx Std. Err. -0.3332 2.1983 > if( mfxLoaded ) { + print( urbinEla( estProbitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 1 ], + xMeanQuad["age"], 1:2, iPos = 0, + model = "lpm", estProbitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 2 ] ) ) + } semEla stdEr -0.3332 1.5544 Warning message: In urbinEla(allCoef = estProbitQuadMfx$mfxest[c("age", "I(age^2)"), 1], allXVal = xMeanQuad["age"], xPos = 1:2, model = "lpm", allCoefVcov = estProbitQuadMfx$mfxest[c("age", "I(age^2)"), 2], iPos = 0) : the returned standard error is likely largely upward biased and, thus, in most cases meaningless; you can provide the full covariance matrix via argument 'allCoefVcov' to avoid this bias or use argument 'xMeanSd' to substantially reduce this bias > if( mfxLoaded ) { + print( urbinEla( estProbitQuadMfx$mfxest[ , 1 ], + xMeanQuad[-1], 2:3, iPos = 0, + model = "lpm", estProbitQuadMfx$mfxest[ , 2 ] ) ) + } semEla stdEr -0.3332 1.5544 Warning message: In urbinEla(allCoef = estProbitQuadMfx$mfxest[, 1], allXVal = xMeanQuad[-1], xPos = 2:3, model = "lpm", allCoefVcov = estProbitQuadMfx$mfxest[, 2], iPos = 0) : the returned standard error is likely largely upward biased and, thus, in most cases meaningless; you can provide the full covariance matrix via argument 'allCoefVcov' to avoid this bias or use argument 'xMeanSd' to substantially reduce this bias > if( mfxLoaded ) { + print( urbinEla( estProbitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 1 ], + xMeanQuad["age"], 1:2, iPos = 0, + model = "lpm", estProbitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 2 ], + xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) ) ) + } semEla stdEr -0.3332 0.1411 > if( mfxLoaded ) { + print( urbinEla( estProbitQuadMfx$mfxest[ , 1 ], + xMeanQuad[-1], 2:3, iPos = 0, + model = "lpm", estProbitQuadMfx$mfxest[ , 2 ], + xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) ) ) + } semEla stdEr -0.3332 0.1411 > > > ### age is interval-coded (age is in the range 30-60) > # create dummy variables for age intervals > Mroz87$age30.37 <- Mroz87$age >= 30 & Mroz87$age <= 37 > Mroz87$age38.44 <- Mroz87$age >= 38 & Mroz87$age <= 44 > Mroz87$age45.52 <- Mroz87$age >= 45 & Mroz87$age <= 52 > Mroz87$age53.60 <- Mroz87$age >= 53 & Mroz87$age <= 60 > all.equal( + Mroz87$age30.37 + Mroz87$age38.44 + Mroz87$age45.52 + Mroz87$age53.60, + rep( 1, nrow( Mroz87 ) ) ) [1] TRUE > # estimation > estProbitInt <- glm( lfp ~ kids + age30.37 + age38.44 + age53.60 + educ, + family = binomial(link = "probit"), + data = Mroz87 ) > summary( estProbitInt ) Call: glm(formula = lfp ~ kids + age30.37 + age38.44 + age53.60 + educ, family = binomial(link = "probit"), data = Mroz87) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.8916 0.2769 -3.22 0.0013 ** kids -0.3775 0.1267 -2.98 0.0029 ** age30.37TRUE 0.1351 0.1270 1.06 0.2876 age38.44TRUE 0.1215 0.1360 0.89 0.3720 age53.60TRUE -0.5142 0.1630 -3.15 0.0016 ** educ 0.1080 0.0214 5.06 4.2e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1029.75 on 752 degrees of freedom Residual deviance: 986.78 on 747 degrees of freedom AIC: 998.8 Number of Fisher Scoring iterations: 4 > # mean values of the explanatory variables > xMeanInt <- c( xMeanLin[1:2], mean( Mroz87$age30.37 ), + mean( Mroz87$age38.44 ), mean( Mroz87$age53.60 ), xMeanLin[4] ) > # semi-elasticity of age without standard errors > urbinElaInt( coef( estProbitInt ), xMeanInt, + c( 3, 4, 0, 5 ), c( 30, 37.5, 44.5, 52.5, 60 ), model = "probit" ) semEla stdEr -0.3754 NA > # semi-elasticities based on numerical derivation > Mroz87Lower <- Mroz87 > Mroz87Lower$age <- Mroz87$age * 0.95 > Mroz87Lower$age30.37 <- Mroz87Lower$age <= 37.5 > Mroz87Lower$age38.44 <- Mroz87Lower$age > 37.5 & Mroz87Lower$age <= 44.5 > Mroz87Lower$age45.52 <- Mroz87Lower$age > 44.5 & Mroz87Lower$age <= 52.5 > Mroz87Lower$age53.60 <- Mroz87Lower$age > 52.5 > all.equal( + Mroz87Lower$age30.37 + Mroz87Lower$age38.44 + Mroz87Lower$age45.52 + + Mroz87Lower$age53.60, rep( 1, nrow( Mroz87 ) ) ) [1] TRUE > Mroz87Upper <- Mroz87 > Mroz87Upper$age <- Mroz87$age * 1.05 > Mroz87Upper$age30.37 <- Mroz87Upper$age <= 37.5 > Mroz87Upper$age38.44 <- Mroz87Upper$age > 37.5 & Mroz87Upper$age <= 44.5 > Mroz87Upper$age45.52 <- Mroz87Upper$age > 44.5 & Mroz87Upper$age <= 52.5 > Mroz87Upper$age53.60 <- Mroz87Upper$age > 52.5 > all.equal( + Mroz87Upper$age30.37 + Mroz87Upper$age38.44 + Mroz87Upper$age45.52 + + Mroz87Upper$age53.60, rep( 1, nrow( Mroz87 ) ) ) [1] TRUE > 10 * mean( predict( estProbitInt, newdata = Mroz87Upper, type = "response" ) - + predict( estProbitInt, newdata = Mroz87Lower, type = "response" ) ) [1] -0.3517 > Mroz87LowerMean <- Mroz87Lower > Mroz87UpperMean <- Mroz87Upper > Mroz87LowerMean$kids <- Mroz87UpperMean$kids <- xMeanInt[ "kids" ] > Mroz87LowerMean$educ <- Mroz87UpperMean$educ <- xMeanInt[ "educ" ] > 10 * mean( predict( estProbitInt, newdata = Mroz87UpperMean, type = "response" ) - + predict( estProbitInt, newdata = Mroz87LowerMean, type = "response" ) ) [1] -0.3664 > # partial derivatives of the semi-elasticity wrt the coefficients > urbinElaInt( coef( estProbitInt ), xMeanInt, + c( 3, 4, 0, 5 ), c( 30, 37.5, 44.5, 52.5, 60 ), model = "probit" )$derivCoef [1] -0.004223 -0.002939 -0.555180 -0.049412 0.545345 -0.051890 > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients > if( maxLikLoaded ) { + print( numericGradient( function( x, ... ){ urbinElaInt( x, ... )$semEla }, + t0 = coef( estProbitInt ), allXVal = xMeanInt, + xPos = c( 3, 4, 0, 5 ), xBound = c( 30, 37.5, 44.5, 52.5, 60 ), + model = "probit" ) ) + } (Intercept) kids age30.37TRUE age38.44TRUE age53.60TRUE educ [1,] -0.004223 -0.002939 -0.5552 -0.04941 0.5453 -0.05189 > # semi-elasticity of age with standard errors (full covariance matrix) > urbinElaInt( coef( estProbitInt ), xMeanInt, + c( 3, 4, 0, 5 ), c( 30, 37.5, 44.5, 52.5, 60 ), + model = "probit", vcov( estProbitInt ) ) semEla stdEr -0.3754 0.1010 > # semi-elasticity of age with standard errors (only standard errors) > urbinElaInt( coef( estProbitInt ), xMeanInt, + c( 3, 4, 0, 5 ), c( 30, 37.5, 44.5, 52.5, 60 ), + model = "probit", sqrt( diag( vcov( estProbitInt ) ) ) ) semEla stdEr -0.3754 0.1137 > # semi-elasticity of age based on partial derivatives calculated by the mfx package > if( mfxLoaded ) { + estProbitIntMfx <- probitmfx( lfp ~ kids + age30.37 + age38.44 + age53.60 + educ, + data = Mroz87 ) + print( urbinElaInt( estProbitIntMfx$mfxest[ 2:4, 1 ], xMeanInt[ 3:5 ], + c( 1, 2, 0, 3 ), iPos = 0, c( 30, 37.5, 44.5, 52.5, 60 ), model = "lpm", + estProbitIntMfx$mfxest[ 2:4, 2 ] ) ) + } semEla stdEr -0.3770 0.1162 > if( mfxLoaded ) { + print( urbinElaInt( estProbitIntMfx$mfxest[ , 1 ], xMeanInt[ -1 ], + c( 2, 3, 0, 4 ), iPos = 0, c( 30, 37.5, 44.5, 52.5, 60 ), model = "lpm", + estProbitIntMfx$mfxest[ , 2 ] ) ) + } semEla stdEr -0.3770 0.1162 > > > ### effect of age changing between discrete intervals > ### if age is used as linear explanatory variable > # mean values of the 'other' explanatory variables > xMeanLinInt <- c( xMeanLin[ 1:2 ], NA, xMeanLin[4] ) > # effects of age changing from the 30-40 interval to the 50-60 interval > # without standard errors > urbinEffInt( coef( estProbitLin ), xMeanLinInt, 3, + c( 30, 40 ), c( 50, 60 ), model = "probit" ) effect stdEr -0.1506 NA > # effects of age changing from the 30-40 interval to the 50-60 interval > # based on predicted values > predict( estProbitLin, + newdata = as.data.frame( t( replace( xMeanLin, 3, 55 ) ) ), + type = "response" ) - + predict( estProbitLin, + newdata = as.data.frame( t( replace( xMeanLin, 3, 35 ) ) ), + type = "response" ) 1 -0.1506 > # partial derivatives of the semi-elasticity wrt the coefficients > urbinEffInt( coef( estProbitLin ), xMeanLinInt, 3, + c( 30, 40 ), c( 50, 60 ), model = "probit" )$derivCoef [1] 0.01969 0.01370 8.65401 0.24193 > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients > if( maxLikLoaded ) { + print( numericGradient( function( x, ... ){ urbinEffInt( x, ... )$effect }, + t0 = coef( estProbitLin ), + allXVal = xMeanLinInt, xPos = 3, + refBound = c( 30, 40 ), intBound = c( 50, 60 ), model = "probit" ) ) + } (Intercept) kids age educ [1,] 0.01969 0.0137 8.654 0.2419 > # effects of age changing from the 30-40 interval to the 50-60 interval > # (full covariance matrix) > urbinEffInt( coef( estProbitLin ), xMeanLinInt, 3, + c( 30, 40 ), c( 50, 60 ), model = "probit", + allCoefVcov = vcov( estProbitLin ) ) effect stdEr -0.15062 0.05457 > # effects of age changing from the 30-40 interval to the 50-60 interval > # (only standard errors) > urbinEffInt( coef( estProbitLin ), xMeanLinInt, 3, + c( 30, 40 ), c( 50, 60 ), model = "probit", + allCoefVcov = sqrt( diag( vcov( estProbitLin ) ) ) ) effect stdEr -0.15062 0.06132 > # semi-elasticity of age based on partial derivative calculated by the mfx package > if( mfxLoaded ) { + print( urbinEffInt( estProbitLinMfx$mfxest[ "age", 1 ], NULL, 1, iPos = 0, + c( 30, 40 ), c( 50, 60 ), model = "lpm", + estProbitLinMfx$mfxest[ "age", 2 ] ) ) + } effect stdEr -0.15041 0.05482 > if( mfxLoaded ) { + print( urbinEffInt( estProbitLinMfx$mfxest[ , 1 ], NULL, 2, iPos = 0, + c( 30, 40 ), c( 50, 60 ), model = "lpm", + estProbitLinMfx$mfxest[ , 2 ] ) ) + } effect stdEr -0.15041 0.05482 > > > ### effect of age changing between discrete intervals > ### if age is used as linear and quadratic explanatory variable > # mean values of the 'other' explanatory variables > xMeanQuadInt <- c( xMeanLin[ 1:2 ], NA, NA, xMeanLin[4] ) > # effects of age changing from the 30-40 interval to the 50-60 interval > # without standard errors > urbinEffInt( coef( estProbitQuad ), xMeanQuadInt, c( 3, 4 ), + c( 30, 40 ), c( 50, 60 ), model = "probit" ) effect stdEr -0.2528 NA > # effects of age changing from the 30-40 interval to the 50-60 interval > # based on predicted values > predict( estProbitQuad, + newdata = as.data.frame( t( replace( xMeanQuad, 3:4, c( 55, 55^2 ) ) ) ), + type = "response" ) - + predict( estProbitQuad, + newdata = as.data.frame( t( replace( xMeanQuad, 3:4, c( 35, 35^2 ) ) ) ), + type = "response" ) 1 -0.2527 > # partial derivatives of the effect wrt the coefficients > urbinEffInt( coef( estProbitQuad ), xMeanQuadInt, c( 3, 4 ), + c( 30, 40 ), c( 50, 60 ), model = "probit" )$derivCoef [1] 2.238e-03 1.557e-03 7.675e+00 6.865e+02 2.750e-02 > # numerically computed partial derivatives of the effect wrt the coefficients > if( maxLikLoaded ) { + print( numericGradient( function( x, ... ){ urbinEffInt( x, ... )$effect }, + t0 = coef( estProbitQuad ), + allXVal = xMeanQuadInt, xPos = c( 3, 4 ), + refBound = c( 30, 40 ), intBound = c( 50, 60 ), model = "probit" ) ) + } (Intercept) kids age I(age^2) educ [1,] 0.002238 0.001557 7.675 686.5 0.0275 > # effects of age changing from the 30-40 interval to the 50-60 interval > # (full covariance matrix) > urbinEffInt( coef( estProbitQuad ), xMeanQuadInt, c( 3, 4 ), + c( 30, 40 ), c( 50, 60 ), model = "probit", + allCoefVcov = vcov( estProbitQuad ) ) effect stdEr -0.25282 0.06148 > # effects of age changing from the 30-40 interval to the 50-60 interval > # (only standard errors) > urbinEffInt( coef( estProbitQuad ), xMeanQuadInt, c( 3, 4 ), + c( 30, 40 ), c( 50, 60 ), model = "probit", + allCoefVcov = sqrt( diag( vcov( estProbitQuad ) ) ) ) effect stdEr -0.2528 0.7332 Warning message: In urbinEffInt(allCoef = coef(estProbitQuad), allXVal = xMeanQuadInt, xPos = c(3, 4), refBound = c(30, 40), intBound = c(50, 60), model = "probit", allCoefVcov = sqrt(diag(vcov(estProbitQuad)))) : the returned standard error is likely largely upward biased and, thus, in most cases meaningless; you can provide the full covariance matrix via argument 'allCoefVcov' to avoid this bias or use argument 'xMeanSd' to substantially reduce this bias > # effects of age changing from the 30-40 interval to the 50-60 interval > # (standard errors + mean value and standard deviation of age) > urbinEffInt( coef( estProbitQuad ), xMeanQuadInt, c( 3, 4 ), + c( 30, 40 ), c( 50, 60 ), model = "probit", + allCoefVcov = sqrt( diag( vcov( estProbitQuad ) ) ), + xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) ) effect stdEr -0.25282 0.07137 > # semi-elasticity of age based on partial derivative calculated by the mfx package > if( mfxLoaded ) { + print( urbinEffInt( estProbitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 1 ], NULL, 1:2, + iPos = 0, c( 30, 40 ), c( 50, 60 ), model = "lpm", + estProbitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 2 ] ) ) + } effect stdEr -0.2530 0.7523 Warning message: In urbinEffInt(allCoef = estProbitQuadMfx$mfxest[c("age", "I(age^2)"), 1], allXVal = NULL, xPos = 1:2, refBound = c(30, 40), intBound = c(50, 60), model = "lpm", allCoefVcov = estProbitQuadMfx$mfxest[c("age", "I(age^2)"), 2], iPos = 0) : the returned standard error is likely largely upward biased and, thus, in most cases meaningless; you can provide the full covariance matrix via argument 'allCoefVcov' to avoid this bias or use argument 'xMeanSd' to substantially reduce this bias > if( mfxLoaded ) { + print( urbinEffInt( estProbitQuadMfx$mfxest[ , 1 ], NULL, 2:3, + iPos = 0, c( 30, 40 ), c( 50, 60 ), model = "lpm", + estProbitQuadMfx$mfxest[ , 2 ] ) ) + } effect stdEr -0.2530 0.7523 Warning message: In urbinEffInt(allCoef = estProbitQuadMfx$mfxest[, 1], allXVal = NULL, xPos = 2:3, refBound = c(30, 40), intBound = c(50, 60), model = "lpm", allCoefVcov = estProbitQuadMfx$mfxest[, 2], iPos = 0) : the returned standard error is likely largely upward biased and, thus, in most cases meaningless; you can provide the full covariance matrix via argument 'allCoefVcov' to avoid this bias or use argument 'xMeanSd' to substantially reduce this bias > if( mfxLoaded ) { + print( urbinEffInt( estProbitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 1 ], NULL, 1:2, + iPos = 0, c( 30, 40 ), c( 50, 60 ), model = "lpm", + estProbitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 2 ], + xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) ) ) + } effect stdEr -0.25303 0.07423 > if( mfxLoaded ) { + print( urbinEffInt( estProbitQuadMfx$mfxest[ , 1 ], NULL, 2:3, + iPos = 0, c( 30, 40 ), c( 50, 60 ), model = "lpm", + estProbitQuadMfx$mfxest[ , 2 ], + xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) ) ) + } effect stdEr -0.25303 0.07423 > > > ### grouping and re-basing categorical variables > ### effects of age changing from the 30-44 category to the 53-60 category > # without standard errors > urbinEffCat( coef( estProbitInt ), xMeanInt, c( 3:5 ), c( -1, -1, 1, 0 ), + model = "probit" ) effect stdEr -0.2524 NA > # effects calculated based on predicted values > names( xMeanInt ) <- sub( "TRUE", "", names( coef( estProbitInt ) ) ) > df30.37 <- df38.44 <- df45.52 <- df53.60 <- as.data.frame( t( xMeanInt ) ) > df30.37[ , 3:5 ] <- c( TRUE, FALSE, FALSE ) > df38.44[ , 3:5 ] <- c( FALSE, TRUE, FALSE ) > df45.52[ , 3:5 ] <- c( FALSE, FALSE, FALSE ) > df53.60[ , 3:5 ] <- c( FALSE, FALSE, TRUE ) > predict( estProbitInt, newdata = df53.60, type = "response" ) - + sum( Mroz87$age30.37 ) / sum( Mroz87$age30.37 + Mroz87$age38.44 ) * + predict( estProbitInt, newdata = df30.37, type = "response" ) - + sum( Mroz87$age38.44 ) / sum( Mroz87$age30.37 + Mroz87$age38.44 ) * + predict( estProbitInt, newdata = df38.44, type = "response" ) 1 -0.2524 > # partial derivatives of the effect wrt the coefficients > urbinEffCat( coef( estProbitInt ), xMeanInt, c( 3:5 ), c( -1, -1, 1, 0 ), + model = "probit" )$derivCoef [1] -0.004644 -0.003231 -0.224001 -0.157068 0.376426 -0.057056 > # numerically computed partial derivatives of the effect wrt the coefficients > if( maxLikLoaded ) { + print( numericGradient( function( x, ... ){ urbinEffCat( x, ... )$effect }, + t0 = coef( estProbitInt ), + allXVal = xMeanInt, xPos = c( 3:5 ), xGroups = c( -1, -1, 1, 0 ), + model = "probit" ) ) + } (Intercept) kids age30.37TRUE age38.44TRUE age53.60TRUE educ [1,] -0.004644 -0.003231 -0.224 -0.1571 0.3764 -0.05706 > # with full covariance matrix > urbinEffCat( coef( estProbitInt ), xMeanInt, c( 3:5 ), c( -1, -1, 1, 0 ), + model = "probit", allCoefVcov = vcov( estProbitInt ) ) effect stdEr -0.25237 0.06538 > # with standard errors only > urbinEffCat( coef( estProbitInt ), xMeanInt, c( 3:5 ), c( -1, -1, 1, 0 ), + model = "probit", allCoefVcov = sqrt( diag( vcov( estProbitInt ) ) ) ) effect stdEr -0.25237 0.07097 > # semi-elasticity of age based on partial derivative calculated by the mfx package > if( mfxLoaded ) { + print( urbinEffCat( estProbitIntMfx$mfxest[ 2:4, 1 ], + xMeanInt[ 3:5 ], c(1:3), iPos = 0, c( -1, -1, 1, 0 ), model = "lpm", + estProbitIntMfx$mfxest[ 2:4, 2 ] ) ) + } effect stdEr -0.25340 0.07249 > if( mfxLoaded ) { + print( urbinEffCat( estProbitIntMfx$mfxest[ , 1 ], + xMeanInt[ -1 ], c(2:4), iPos = 0, c( -1, -1, 1, 0 ), model = "lpm", + estProbitIntMfx$mfxest[ , 2 ] ) ) + } effect stdEr -0.25340 0.07249 > > > ### effects of age changing from the 53-60 category to the 38-52 category > # without standard errors > urbinEffCat( coef( estProbitInt ), xMeanInt, c( 3:5 ), c( 0, 1, -1, 1 ), + model = "probit" ) effect stdEr 0.2228 NA > # effects calculated based on predicted values > sum( Mroz87$age38.44 ) / sum( Mroz87$age38.44 + Mroz87$age45.52 ) * + predict( estProbitInt, newdata = df38.44, type = "response" ) + + sum( Mroz87$age45.52 ) / sum( Mroz87$age38.44 + Mroz87$age45.52 ) * + predict( estProbitInt, newdata = df45.52, type = "response" ) - + predict( estProbitInt, newdata = df53.60, type = "response" ) 1 0.2227 > # partial derivatives of the effect wrt the coefficients > urbinEffCat( coef( estProbitInt ), xMeanInt, c( 3:5 ), c( 0, 1, -1, 1 ), + model = "probit" )$derivCoef [1] 0.012449 0.008663 0.000000 0.168992 -0.376426 0.152958 > # numerically computed partial derivatives of the effect wrt the coefficients > if( maxLikLoaded ) { + print( numericGradient( function( x, ... ){ urbinEffCat( x, ... )$effect }, + t0 = coef( estProbitInt ), + allXVal = xMeanInt, xPos = c( 3:5 ), xGroups = c( 0, 1, -1, 1 ), + model = "probit" ) ) + } (Intercept) kids age30.37TRUE age38.44TRUE age53.60TRUE educ [1,] 0.01245 0.008663 0 0.169 -0.3764 0.153 > # with full covariance matrix > urbinEffCat( coef( estProbitInt ), xMeanInt, c( 3:5 ), c( 0, 1, -1, 1 ), + model = "probit", allCoefVcov = vcov( estProbitInt ) ) effect stdEr 0.2228 0.0607 > # with standard errors only > urbinEffCat( coef( estProbitInt ), xMeanInt, c( 3:5 ), c( 0, 1, -1, 1 ), + model = "probit", allCoefVcov = sqrt( diag( vcov( estProbitInt ) ) ) ) effect stdEr 0.22284 0.06572 > # semi-elasticity of age based on partial derivative calculated by the mfx package > if( mfxLoaded ) { + print( urbinEffCat( estProbitIntMfx$mfxest[ 2:4, 1 ], xMeanInt[ 3:5 ], c(1:3), + c( 0, 1, -1, 1 ), iPos = 0, model = "lpm", estProbitIntMfx$mfxest[ 2:4, 2 ] ) ) + } effect stdEr 0.22346 0.06683 > if( mfxLoaded ) { + print( urbinEffCat( estProbitIntMfx$mfxest[ , 1 ], xMeanInt[ -1 ], c(2:4), + c( 0, 1, -1, 1 ), iPos = 0, model = "lpm", estProbitIntMfx$mfxest[ , 2 ] ) ) + } effect stdEr 0.22346 0.06683 > > > proc.time() user system elapsed 1.81 0.18 1.98