Running 'logit.R' [160s]
  Running 'logitEffCat.R' [0s]
  Comparing 'logitEffCat.Rout' to 'logitEffCat.Rout.save' ... OK
  Running 'logitEffInt.R' [0s]
  Comparing 'logitEffInt.Rout' to 'logitEffInt.Rout.save' ... OK
  Running 'logitEla.R' [0s]
  Comparing 'logitEla.Rout' to 'logitEla.Rout.save' ... OK
  Running 'logitElaInt.R' [174s]
  Running 'lpm.R' [163s]
  Running 'mlogit.R' [5s]
  Comparing 'mlogit.Rout' to 'mlogit.Rout.save' ... OK
  Running 'mvprobit.R' [18s]
  Comparing 'mvprobit.Rout' to 'mvprobit.Rout.save' ... OK
  Running 'oprobit.R' [2s]
  Comparing 'oprobit.Rout' to 'oprobit.Rout.save' ... OK
  Running 'probit.R' [162s]
Running the tests in 'tests/logit.R' failed.
Complete output:
  > 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
  > estLogitLin <- glm( lfp ~ kids + age + educ,
  +   family = binomial(link = "logit"),
  +   data = Mroz87 )
  > summary( estLogitLin )
  
  Call:
  glm(formula = lfp ~ kids + age + educ, family = binomial(link = "logit"), 
      data = Mroz87)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept)  -0.1144     0.7346   -0.16   0.8763    
  kids         -0.5035     0.1993   -2.53   0.0115 *  
  age          -0.0311     0.0114   -2.73   0.0063 ** 
  educ          0.1690     0.0351    4.82  1.4e-06 ***
  ---
  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.53  on 749  degrees of freedom
  AIC: 1002
  
  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( estLogitLin ), xMeanLin, xPos = 3, model = "logit" )
   semEla   stdEr 
  -0.3237      NA 
  > # semi-elasticity of age based on numerical derivation
  > 100 * ( predict( estLogitLin,
  +   newdata = as.data.frame( t( xMeanLin * c( 1, 1, 1.005, 1 ) ) ),
  +   type = "response" ) -
  +     predict( estLogitLin,
  +       newdata = as.data.frame( t( xMeanLin * c( 1, 1, 0.995, 1 ) ) ),
  +       type = "response" ) )
        1 
  -0.3237 
  > # partial derivatives of the semi-elasticity wrt the coefficients
  > urbinEla( coef( estLogitLin ), xMeanLin, 3, seSimplify = FALSE,
  +   model = "logit" )$derivCoef
               kids      age     educ 
   0.04658  0.03242 12.39575  0.57236 
  > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEla( x, ... )$semEla },
  +     t0 = coef( estLogitLin ),
  +     allXVal = xMeanLin, xPos = 3, model = "logit" ) )
  + }
       (Intercept)    kids  age   educ
  [1,]     0.04658 0.03242 12.4 0.5724
  > # simplified partial derivatives of the semi-elasticity wrt the coefficients
  > urbinEla( coef( estLogitLin ), xMeanLin, 3, model = "logit",
  +   seSimplify = TRUE )$derivCoef
  [1]  0.00  0.00 10.41  0.00
  > # semi-elasticity of age with standard errors (full covariance matrix)
  > urbinEla( coef( estLogitLin ), xMeanLin, 3, model = "logit",
  +   vcov( estLogitLin ) )
   semEla   stdEr 
  -0.3237  0.1183 
  > # semi-elasticity of age with standard errors (only standard errors)
  > urbinEla( coef( estLogitLin ), xMeanLin, 3, model = "logit",
  +   sqrt( diag( vcov( estLogitLin ) ) ), seSimplify = FALSE )
   semEla   stdEr 
  -0.3237  0.1465 
  Warning message:
  In urbinEla(coef(estLogitLin), xMeanLin, 3, model = "logit", sqrt(diag(vcov(estLogitLin))),  :
    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( estLogitLin ), xMeanLin, 3, model = "logit",
  +   sqrt( diag( vcov( estLogitLin ) ) ) )
   semEla   stdEr 
  -0.3237  0.1184 
  > # semi-elasticity of age based on partial derivative calculated by the mfx package
  > if( mfxLoaded ) {
  +   estLogitLinMfx <- logitmfx( lfp ~ kids + age + educ, data = Mroz87 )
  +   print( estLogitLinMfx$mfxest[ "age", 1:2 ] * xMeanLin[ "age" ] )
  + }
      dF/dx Std. Err. 
    -0.3237    0.1183 
  > if( mfxLoaded ) {
  +   print( urbinEla( estLogitLinMfx$mfxest[ "age", 1 ], xMeanLin["age"], 1, iPos = 0,
  +     model = "lpm", estLogitLinMfx$mfxest[ "age", 2 ] ) )
  + }
   semEla   stdEr 
  -0.3237  0.1183 
  > if( mfxLoaded ) {
  +   print( urbinEla( estLogitLinMfx$mfxest[ , 1 ], xMeanLin[-1], 2, iPos = 0,
  +     model = "lpm", estLogitLinMfx$mfxest[ , 2 ] ) )
  + }
   semEla   stdEr 
  -0.3237  0.1183 
  > 
  > 
  > ### quadratic in age
  > estLogitQuad <- glm( lfp ~ kids + age + I(age^2) + educ,
  +   family = binomial(link = "logit"),
  +   data = Mroz87 )
  > summary( estLogitQuad )
  
  Call:
  glm(formula = lfp ~ kids + age + I(age^2) + educ, family = binomial(link = "logit"), 
      data = Mroz87)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept) -6.94366    2.27412   -3.05  0.00226 ** 
  kids        -0.70271    0.21334   -3.29  0.00099 ***
  age          0.30802    0.10737    2.87  0.00412 ** 
  I(age^2)    -0.00400    0.00126   -3.17  0.00153 ** 
  educ         0.17229    0.03545    4.86  1.2e-06 ***
  ---
  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:  983.26  on 748  degrees of freedom
  AIC: 993.3
  
  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( estLogitQuad ), xMeanQuad, c( 3, 4 ), model = "logit" )
   semEla   stdEr 
  -0.3172      NA 
  > # semi-elasticity of age based on numerical derivation
  > 100 * ( predict( estLogitQuad,
  +   newdata = as.data.frame( t( xMeanQuad * c( 1, 1, 1.005, 1.005^2, 1 ) ) ),
  +   type = "response" ) -
  +     predict( estLogitQuad,
  +       newdata = as.data.frame( t( xMeanQuad * c( 1, 1, 0.995, 0.995^2, 1 ) ) ),
  +       type = "response" ) )
        1 
  -0.3172 
  > # partial derivatives of the semi-elasticity wrt the coefficients
  > urbinEla( coef( estLogitQuad ), xMeanQuad, c( 3, 4 ), model = "logit",
  +   seSimplify = FALSE )$derivCoef
                 kids       age       age      educ 
    0.08535   0.05939  13.49480 993.64782   1.04864 
  > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEla( x, ... )$semEla },
  +     t0 = coef( estLogitQuad ),
  +     allXVal = xMeanQuad, xPos = c( 3, 4 ), model = "logit" ) )
  + }
       (Intercept)    kids   age I(age^2)  educ
  [1,]     0.08535 0.05939 13.49    993.6 1.049
  > # simplified partial derivatives of the semi-elasticity wrt the coefficients
  > urbinEla( coef( estLogitQuad ), xMeanQuad, c( 3, 4 ), model = "logit",
  +   seSimplify = TRUE )$derivCoef
  [1]   0.000   0.000   9.864 839.216   0.000
  > # semi-elasticity of age with standard errors (full covariance matrix)
  > urbinEla( coef( estLogitQuad ), xMeanQuad, c( 3, 4 ), model = "logit",
  +   vcov( estLogitQuad ) )
   semEla   stdEr 
  -0.3172  0.1125 
  > # semi-elasticity of age with standard errors (only standard errors)
  > urbinEla( coef( estLogitQuad ), xMeanQuad, c( 3, 4 ), model = "logit",
  +   sqrt( diag( vcov( estLogitQuad ) ) ), seSimplify = FALSE )
   semEla   stdEr 
  -0.3172  1.9265 
  Warning messages:
  1: In urbinEla(coef(estLogitQuad), xMeanQuad, c(3, 4), model = "logit",  :
    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(estLogitQuad), allXVal = xMeanQuad, xPos = c(3, 4), model = "logit", allCoefVcov = sqrt(diag(vcov(estLogitQuad))), 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( estLogitQuad ), xMeanQuad, c( 3, 4 ), model = "logit",
  +   sqrt( diag( vcov( estLogitQuad ) ) ) )
   semEla   stdEr 
  -0.3172  1.4979 
  Warning message:
  In urbinEla(allCoef = coef(estLogitQuad), allXVal = xMeanQuad, xPos = c(3, 4), model = "logit", allCoefVcov = sqrt(diag(vcov(estLogitQuad)))) :
    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( estLogitQuad ), xMeanQuad, c( 3, 4 ), model = "logit",
  +   sqrt( diag( vcov( estLogitQuad ) ) ),
  +   xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ),
  +   seSimplify = FALSE )
    semEla    stdEr 
  -0.31715  0.09913 
  Warning message:
  In urbinEla(coef(estLogitQuad), xMeanQuad, c(3, 4), model = "logit",  :
    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( estLogitQuad ), xMeanQuad, c( 3, 4 ), model = "logit",
  +   sqrt( diag( vcov( estLogitQuad ) ) ),
  +   xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) )
   semEla   stdEr 
  -0.3172  0.1360 
  > # 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 ) {
  +   estLogitQuadMfx <- logitmfx( lfp ~ kids + age + I(age^2) + educ, data = Mroz87 )
  +   print( estLogitQuadMfx$mfxest[ "age", 1:2 ] * xMeanQuad[ "age" ] +
  +     2 * estLogitQuadMfx$mfxest[ "I(age^2)", 1:2 ] * xMeanQuad[ "age" ]^2 )
  + }
      dF/dx Std. Err. 
    -0.3347    2.2359 
  > if( mfxLoaded ) {
  +   print( urbinEla( estLogitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 1 ],
  +     xMeanQuad["age"], 1:2, iPos = 0,
  +     model = "lpm", estLogitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 2 ] ) )
  + }
   semEla   stdEr 
  -0.3347  1.5810 
  Warning message:
  In urbinEla(allCoef = estLogitQuadMfx$mfxest[c("age", "I(age^2)"), 1], allXVal = xMeanQuad["age"], xPos = 1:2, model = "lpm", allCoefVcov = estLogitQuadMfx$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( estLogitQuadMfx$mfxest[  , 1 ],
  +     xMeanQuad[-1], 2:3, iPos = 0,
  +     model = "lpm", estLogitQuadMfx$mfxest[ , 2 ] ) )
  + }
   semEla   stdEr 
  -0.3347  1.5810 
  Warning message:
  In urbinEla(allCoef = estLogitQuadMfx$mfxest[, 1], allXVal = xMeanQuad[-1], xPos = 2:3, model = "lpm", allCoefVcov = estLogitQuadMfx$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( estLogitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 1 ],
  +     xMeanQuad["age"], 1:2, iPos = 0,
  +     model = "lpm", estLogitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 2 ],
  +     xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) ) )
  + }
   semEla   stdEr 
  -0.3347  0.1435 
  > if( mfxLoaded ) {
  +   print( urbinEla( estLogitQuadMfx$mfxest[  , 1 ],
  +     xMeanQuad[-1], 2:3, iPos = 0,
  +     model = "lpm", estLogitQuadMfx$mfxest[ , 2 ],
  +     xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) ) )
  + }
   semEla   stdEr 
  -0.3347  0.1435 
  > 
  > 
  > ### 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
  > estLogitInt <- glm( lfp ~ kids + age30.37 + age38.44 + age53.60 + educ,
  +   family = binomial(link = "logit"),
  +   data = Mroz87 )
  > summary( estLogitInt )
  
  Call:
  glm(formula = lfp ~ kids + age30.37 + age38.44 + age53.60 + educ, 
      family = binomial(link = "logit"), data = Mroz87)
  
  Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
  (Intercept)   -1.4441     0.4548   -3.18   0.0015 ** 
  kids          -0.6045     0.2080   -2.91   0.0037 ** 
  age30.37TRUE   0.2208     0.2060    1.07   0.2839    
  age38.44TRUE   0.1973     0.2203    0.90   0.3704    
  age53.60TRUE  -0.8271     0.2653   -3.12   0.0018 ** 
  educ           0.1743     0.0354    4.93  8.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:  987.04  on 747  degrees of freedom
  AIC: 999
  
  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( estLogitInt ), xMeanInt,
  +   c( 3, 4, 0, 5 ), c( 30, 37.5, 44.5, 52.5, 60 ), model = "logit" )
   semEla   stdEr 
  -0.3777      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( estLogitInt, newdata = Mroz87Upper, type = "response" ) -
  +   predict( estLogitInt, newdata = Mroz87Lower, type = "response" ) )
  [1] -0.351
  > Mroz87LowerMean <- Mroz87Lower
  > Mroz87UpperMean <- Mroz87Upper
  > Mroz87LowerMean$kids <- Mroz87UpperMean$kids <- xMeanInt[ "kids" ]
  > Mroz87LowerMean$educ <- Mroz87UpperMean$educ <- xMeanInt[ "educ" ]
  > 10 * mean( predict( estLogitInt, newdata = Mroz87UpperMean, type = "response" ) -
  +     predict( estLogitInt, newdata = Mroz87LowerMean, type = "response" ) )
  [1] -0.3688
  > # partial derivatives of the semi-elasticity wrt the coefficients
  > urbinElaInt( coef( estLogitInt ), xMeanInt,
  +   c( 3, 4, 0, 5 ), c( 30, 37.5, 44.5, 52.5, 60 ), model = "logit" )$derivCoef
  [1] -0.003229 -0.002247 -0.343126 -0.030585  0.336143 -0.039678
  > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinElaInt( x, ... )$semEla },
  +     t0 = coef( estLogitInt ), allXVal = xMeanInt,
  +     xPos = c( 3, 4, 0, 5 ), xBound = c( 30, 37.5, 44.5, 52.5, 60 ),
  +     model = "logit" ) )
  + }
       (Intercept)      kids age30.37TRUE age38.44TRUE age53.60TRUE     educ
  [1,]   -0.003229 -0.002247      -0.3431     -0.03058       0.3361 -0.03968
  > # semi-elasticity of age with standard errors (full covariance matrix)
  > urbinElaInt( coef( estLogitInt ), xMeanInt,
  +   c( 3, 4, 0, 5 ), c( 30, 37.5, 44.5, 52.5, 60 ), model = "logit",
  +   allCoefVcov = vcov( estLogitInt ) )
   semEla   stdEr 
  -0.3777  0.1017 
  > # semi-elasticity of age with standard errors (only standard errors)
  > urbinElaInt( coef( estLogitInt ), xMeanInt,
  +   c( 3, 4, 0, 5 ), c( 30, 37.5, 44.5, 52.5, 60 ), model = "logit",
  +   allCoefVcov = sqrt( diag( vcov( estLogitInt ) ) ) )
   semEla   stdEr 
  -0.3777  0.1140 
  > # semi-elasticity of age based on partial derivatives calculated by the mfx package
  > if( mfxLoaded ) {
  +   estLogitIntMfx <- logitmfx( lfp ~ kids + age30.37 + age38.44 + age53.60 + educ,
  +     data = Mroz87 )
  +   print( urbinElaInt( estLogitIntMfx$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",
  +     estLogitIntMfx$mfxest[ 2:4, 2 ] ) )
  + }
   semEla   stdEr 
  -0.3799  0.1173 
  > if( mfxLoaded ) {
  +   print( urbinElaInt( estLogitIntMfx$mfxest[ , 1 ], xMeanInt[ -1 ],
  +     c( 2, 3, 0, 4 ), iPos = 0, c( 30, 37.5, 44.5, 52.5, 60 ), model = "lpm",
  +     estLogitIntMfx$mfxest[ , 2 ] ) )
  + }
   semEla   stdEr 
  -0.3799  0.1173 
  > 
  > 
  > ### 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( estLogitLin ), allXVal = xMeanLinInt, xPos = 3,
  +   refBound = c( 30, 40 ), intBound = c( 50, 60 ), model = "logit" )
   effect   stdEr 
  -0.1525      NA 
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # based on predicted values
  > predict( estLogitLin,
  +   newdata = as.data.frame( t( replace( xMeanLin, 3, 55 ) ) ),
  +   type = "response" ) -
  +   predict( estLogitLin,
  +     newdata = as.data.frame( t( replace( xMeanLin, 3, 35 ) ) ),
  +     type = "response" )
        1 
  -0.1525 
  > # partial derivatives of the semi-elasticity wrt the coefficients
  > urbinEffInt( coef( estLogitLin ), xMeanLinInt, 3,
  +   c( 30, 40 ), c( 50, 60 ), model = "logit" )$derivCoef
  [1] 0.01582 0.01101 5.54183 0.19437
  > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffInt( x, ... )$effect },
  +     t0 = coef( estLogitLin ),
  +     allXVal = xMeanLinInt, xPos = 3,
  +     refBound = c( 30, 40 ), intBound = c( 50, 60 ), model = "logit" ) )
  + }
       (Intercept)    kids   age   educ
  [1,]     0.01582 0.01101 5.542 0.1944
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # (full covariance matrix)
  > urbinEffInt( coef( estLogitLin ), xMeanLinInt, 3,
  +   c( 30, 40 ), c( 50, 60 ), model = "logit", allCoefVcov = vcov( estLogitLin ) )
    effect    stdEr 
  -0.15247  0.05533 
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # (only standard errors)
  > urbinEffInt( coef( estLogitLin ), allXVal = xMeanLinInt, xPos = 3,
  +   refBound = c( 30, 40 ), intBound = c( 50, 60 ), model = "logit",
  +   allCoefVcov = sqrt( diag( vcov( estLogitLin ) ) ) )
    effect    stdEr 
  -0.15247  0.06446 
  > # semi-elasticity of age based on partial derivative calculated by the mfx package
  > if( mfxLoaded ) {
  +   print( urbinEffInt( estLogitLinMfx$mfxest[ "age", 1 ], NULL, 1, iPos = 0,
  +     c( 30, 40 ), c( 50, 60 ), model = "lpm",
  +     estLogitLinMfx$mfxest[ "age", 2 ] ) )
  + }
    effect    stdEr 
  -0.15219  0.05564 
  > if( mfxLoaded ) {
  +   print( urbinEffInt( estLogitLinMfx$mfxest[ , 1 ], NULL, 2, iPos = 0,
  +     c( 30, 40 ), c( 50, 60 ), model = "lpm",
  +     estLogitLinMfx$mfxest[ , 2 ] ) )
  + }
    effect    stdEr 
  -0.15219  0.05564 
  > 
  > 
  > ### 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( estLogitQuad ), allXVal = xMeanQuadInt,
  +   xPos = c( 3, 4 ), refBound = c( 30, 40 ), intBound = c( 50, 60 ),
  +   model = "logit" )
   effect   stdEr 
  -0.2535      NA 
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # based on predicted values
  > predict( estLogitQuad,
  +   newdata = as.data.frame( t( replace( xMeanQuad, 3:4, c( 55, 55^2 ) ) ) ),
  +   type = "response" ) -
  +   predict( estLogitQuad,
  +     newdata = as.data.frame( t( replace( xMeanQuad, 3:4, c( 35, 35^2 ) ) ) ),
  +     type = "response" )
        1 
  -0.2534 
  > # partial derivatives of the effect wrt the coefficients
  > urbinEffInt( coef( estLogitQuad ), xMeanQuadInt, c( 3, 4 ),
  +   c( 30, 40 ), c( 50, 60 ), model = "logit" )$derivCoef
  [1] 1.806e-03 1.257e-03 4.760e+00 4.249e+02 2.219e-02
  > # numerically computed partial derivatives of the effect wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffInt( x, ... )$effect },
  +     t0 = coef( estLogitQuad ),
  +     allXVal = xMeanQuadInt, xPos = c( 3, 4 ),
  +     refBound = c( 30, 40 ), intBound = c( 50, 60 ), model = "logit" ) )
  + }
       (Intercept)     kids  age I(age^2)    educ
  [1,]    0.001806 0.001257 4.76    424.9 0.02219
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # (full covariance matrix)
  > urbinEffInt( coef( estLogitQuad ), xMeanQuadInt, c( 3, 4 ),
  +   c( 30, 40 ), c( 50, 60 ), model = "logit",
  +   allCoefVcov = vcov( estLogitQuad ) )
    effect    stdEr 
  -0.25354  0.06224 
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # (only standard errors)
  > urbinEffInt( coef( estLogitQuad ), allXVal = xMeanQuadInt,
  +   xPos = c( 3, 4 ), refBound = c( 30, 40 ), intBound = c( 50, 60 ),
  +   model = "logit", allCoefVcov = sqrt( diag( vcov( estLogitQuad ) ) ) )
   effect   stdEr 
  -0.2535  0.7408 
  Warning message:
  In urbinEffInt(allCoef = coef(estLogitQuad), allXVal = xMeanQuadInt, xPos = c(3, 4), refBound = c(30, 40), intBound = c(50, 60), model = "logit", allCoefVcov = sqrt(diag(vcov(estLogitQuad)))) :
    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( estLogitQuad ), xMeanQuadInt, c( 3, 4 ),
  +   c( 30, 40 ), c( 50, 60 ), model = "logit",
  +   allCoefVcov = sqrt( diag( vcov( estLogitQuad ) ) ),
  +   xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) )
    effect    stdEr 
  -0.25354  0.07205 
  > # semi-elasticity of age based on partial derivative calculated by the mfx package
  > if( mfxLoaded ) {
  +   print( urbinEffInt( estLogitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 1 ], NULL, 1:2,
  +     iPos = 0, c( 30, 40 ), c( 50, 60 ), model = "lpm",
  +     estLogitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 2 ] ) )
  + }
   effect   stdEr 
  -0.2538  0.7652 
  Warning message:
  In urbinEffInt(allCoef = estLogitQuadMfx$mfxest[c("age", "I(age^2)"), 1], allXVal = NULL, xPos = 1:2, refBound = c(30, 40), intBound = c(50, 60), model = "lpm", allCoefVcov = estLogitQuadMfx$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( estLogitQuadMfx$mfxest[ , 1 ], NULL, 2:3,
  +     iPos = 0, c( 30, 40 ), c( 50, 60 ), model = "lpm",
  +     estLogitQuadMfx$mfxest[ , 2 ] ) )
  + }
   effect   stdEr 
  -0.2538  0.7652 
  Warning message:
  In urbinEffInt(allCoef = estLogitQuadMfx$mfxest[, 1], allXVal = NULL, xPos = 2:3, refBound = c(30, 40), intBound = c(50, 60), model = "lpm", allCoefVcov = estLogitQuadMfx$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( estLogitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 1 ], NULL, 1:2,
  +     iPos = 0, c( 30, 40 ), c( 50, 60 ), model = "lpm",
  +     estLogitQuadMfx$mfxest[ c( "age", "I(age^2)" ), 2 ],
  +     xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) ) )
  + }
    effect    stdEr 
  -0.25377  0.07577 
  > if( mfxLoaded ) {
  +   print( urbinEffInt( estLogitQuadMfx$mfxest[ , 1 ], NULL, 2:3,
  +     iPos = 0, c( 30, 40 ), c( 50, 60 ), model = "lpm",
  +     estLogitQuadMfx$mfxest[ , 2 ],
  +     xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) ) )
  + }
    effect    stdEr 
  -0.25377  0.07577 
  > 
  > 
  > ### 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( estLogitInt ), xMeanInt, c( 3:5 ), c( -1, -1, 1, 0 ),
  +   model = "logit" )
   effect   stdEr 
  -0.2538      NA 
  > # effects calculated based on predicted values
  > names( xMeanInt ) <- sub( "TRUE", "", names( coef( estLogitInt ) ) )
  > 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( estLogitInt, newdata = df53.60, type = "response" ) -
  +   sum( Mroz87$age30.37 ) / sum( Mroz87$age30.37 + Mroz87$age38.44 ) *
  +   predict( estLogitInt, newdata = df30.37, type = "response" ) -
  +   sum( Mroz87$age38.44 ) / sum( Mroz87$age30.37 + Mroz87$age38.44 ) *
  +   predict( estLogitInt, newdata = df38.44, type = "response" )
        1 
  -0.2538 
  > # partial derivatives of the effect wrt the coefficients
  > urbinEffCat( coef( estLogitInt ), xMeanInt, c( 3:5 ),
  +   c( -1, -1, 1, 0 ), model = "logit" )$derivCoef
  [1] -0.003642 -0.002534 -0.138529 -0.097136  0.232023 -0.044745
  > # numerically computed partial derivatives of the effect wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffCat( x, ... )$effect },
  +     t0 = coef( estLogitInt ),
  +     allXVal = xMeanInt, xPos = c( 3:5 ), xGroups = c( -1, -1, 1, 0 ),
  +     model = "logit" ) )
  + }
       (Intercept)      kids age30.37TRUE age38.44TRUE age53.60TRUE     educ
  [1,]   -0.003642 -0.002534      -0.1385     -0.09714        0.232 -0.04474
  > # with full covariance matrix
  > urbinEffCat( coef( estLogitInt ), xMeanInt, c( 3:5 ), c( -1, -1, 1, 0 ),
  +   allCoefVcov = vcov( estLogitInt ), model = "logit" )
    effect    stdEr 
  -0.25381  0.06583 
  > # with standard errors only
  > urbinEffCat( coef( estLogitInt ), xMeanInt, c( 3:5 ), c( -1, -1, 1, 0 ),
  +   allCoefVcov = sqrt( diag( vcov( estLogitInt ) ) ), model = "logit" )
    effect    stdEr 
  -0.25381  0.07119 
  > # semi-elasticity of age based on partial derivative calculated by the mfx package
  > if( mfxLoaded ) {
  +   print( urbinEffCat( estLogitIntMfx$mfxest[ 2:4, 1 ], xMeanInt[ 3:5 ], c(1:3),
  +     c( -1, -1, 1, 0 ), iPos = 0, model = "lpm", estLogitIntMfx$mfxest[ 2:4, 2 ] ) )
  + }
    effect    stdEr 
  -0.25514  0.07315 
  > if( mfxLoaded ) {
  +   print( urbinEffCat( estLogitIntMfx$mfxest[ , 1 ], xMeanInt[ -1 ], c(2:4),
  +     c( -1, -1, 1, 0 ), iPos = 0, model = "lpm", estLogitIntMfx$mfxest[ , 2 ] ) )
  + }
    effect    stdEr 
  -0.25514  0.07315 
  > 
  > 
  > ### effects of age changing from the 53-60 category to the 38-52 category
  > # without standard errors
  > urbinEffCat( coef( estLogitInt ), xMeanInt, c( 3:5 ), c( 0, 1, -1, 1 ),
  +   model = "logit" )
  effect  stdEr 
  0.2239     NA 
  > # effects calculated based on predicted values
  > sum( Mroz87$age38.44 ) / sum( Mroz87$age38.44 + Mroz87$age45.52 ) *
  +   predict( estLogitInt, newdata = df38.44, type = "response" ) +
  +   sum( Mroz87$age45.52 ) / sum( Mroz87$age38.44 + Mroz87$age45.52 ) *
  +   predict( estLogitInt, newdata = df45.52, type = "response" ) -
  +   predict( estLogitInt, newdata = df53.60, type = "response" )
       1 
  0.2236 
  > # partial derivatives of the effect wrt the coefficients
  > urbinEffCat( coef( estLogitInt ), xMeanInt, c( 3:5 ),
  +   c( 0, 1, -1, 1 ), model = "logit" )$derivCoef
  [1]  0.009917  0.006901  0.000000  0.105139 -0.232023  0.121846
  > # numerically computed partial derivatives of the effect wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffCat( x, ... )$effect },
  +     t0 = coef( estLogitInt ),
  +     allXVal = xMeanInt, xPos = c( 3:5 ), xGroups = c( 0, 1, -1, 1 ),
  +     model = "logit" ) )
  + }
       (Intercept)     kids age30.37TRUE age38.44TRUE age53.60TRUE   educ
  [1,]    0.009917 0.006901            0       0.1051       -0.232 0.1218
  > # with full covariance matrix
  > urbinEffCat( coef( estLogitInt ), xMeanInt, c( 3:5 ), c( 0, 1, -1, 1 ),
  +   allCoefVcov = vcov( estLogitInt ), model = "logit" )
   effect   stdEr 
  0.22386 0.06109 
  > # with standard errors only
  > urbinEffCat( coef( estLogitInt ), xMeanInt, c( 3:5 ), c( 0, 1, -1, 1 ),
  +   allCoefVcov = sqrt( diag( vcov( estLogitInt ) ) ), model = "logit" )
   effect   stdEr 
  0.22386 0.06609 
  > # semi-elasticity of age based on partial derivative calculated by the mfx package
  > if( mfxLoaded ) {
  +   print( urbinEffCat( estLogitIntMfx$mfxest[ 2:4, 1 ], xMeanInt[ 3:5 ], c(1:3),
  +     c( 0, 1, -1, 1 ), iPos = 0, model = "lpm", estLogitIntMfx$mfxest[ 2:4, 2 ] ) )
  + }
   effect   stdEr 
  0.22466 0.06745 
  > if( mfxLoaded ) {
  +   print( urbinEffCat( estLogitIntMfx$mfxest[ , 1 ], xMeanInt[ -1 ], c(2:4),
  +     c( 0, 1, -1, 1 ), iPos = 0, model = "lpm", estLogitIntMfx$mfxest[ , 2 ] ) )
  + }
   effect   stdEr 
  0.22466 0.06745 
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.51    0.25    2.76 
Running the tests in 'tests/logitElaInt.R' failed.
Complete output:
  > library( "urbin" )
  > 
  > # Example
  > ela8a <- urbinElaInt( 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 ),
  +   model = "logit" )
  > ela8a 
       semEla       stdEr 
  0.008228764          NA 
  > 
  > ela8b <- urbinElaInt( 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 ),
  +   allCoefVcov = c( 0.003, 0.045, 0.007, 0.009 ),
  +   model = "logit" )
  > ela8b 
       semEla       stdEr 
  0.008228764 0.003959305 
  > 
  > # Example
  > ela9a <- urbinElaInt( 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, 
  +   model = "mlogit" )
  > ela9a
      semEla      stdEr 
  0.04186674         NA 
  > 
  > ela9b <- urbinElaInt( 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, 
  +   allCoefVcov = c( 0.003, 0.045, 0.007, 0.009, 0.0008, 0.9 ),
  +   model = "mlogit" )
  > ela9b
      semEla      stdEr 
  0.04186674 0.09570526 
  > 
  > proc.time()
     user  system elapsed 
     0.15    0.06    0.20 
Running the tests in 'tests/lpm.R' failed.
Complete output:
  > 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/
  > 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
  > estLpmLin <- lm( lfp ~ kids + age + educ,
  +   data = Mroz87 )
  > summary( estLpmLin )
  
  Call:
  lm(formula = lfp ~ kids + age + educ, data = Mroz87)
  
  Residuals:
     Min     1Q Median     3Q    Max 
  -0.861 -0.530  0.275  0.431  0.711 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept)  0.47432    0.17069    2.78   0.0056 ** 
  kids        -0.11481    0.04564   -2.52   0.0121 *  
  age         -0.00716    0.00262   -2.73   0.0064 ** 
  educ         0.03893    0.00781    4.98  7.8e-07 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 0.485 on 749 degrees of freedom
  Multiple R-squared:  0.0466,	Adjusted R-squared:  0.0427 
  F-statistic: 12.2 on 3 and 749 DF,  p-value: 8.5e-08
  
  > # mean values of the explanatory variables
  > xMeanLin <- c( 1, colMeans( Mroz87[ , c( "kids", "age", "educ" ) ] ) )
  > # semi-elasticity of age without standard errors
  > urbinEla( coef( estLpmLin )[ "age" ], xMeanLin[ "age" ], xPos = 1, iPos = 0,
  +   model = "lpm" )
   semEla   stdEr 
  -0.3044      NA 
  > urbinEla( coef( estLpmLin ), xMeanLin, xPos = 3, model = "lpm" )
   semEla   stdEr 
  -0.3044      NA 
  > # semi-elasticity of age based on numerical derivation
  > 100 * ( predict( estLpmLin,
  +   newdata = as.data.frame( t( xMeanLin * c( 1, 1, 1.005, 1 ) ) ) ) -
  +     predict( estLpmLin,
  +       newdata = as.data.frame( t( xMeanLin * c( 1, 1, 0.995, 1 ) ) ) ) )
        1 
  -0.3044 
  > # partial derivatives of the semi-elasticity wrt the coefficients
  > urbinEla( coef( estLpmLin )["age"], xMeanLin["age"], xPos = 1, iPos = 0,
  +   model = "lpm" )$derivCoef
  [1] 42.54
  > urbinEla( coef( estLpmLin ), xMeanLin, xPos = 3, model = "lpm" )$derivCoef
  [1]  0.00  0.00 42.54  0.00
  > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEla( x, ... )$semEla },
  +     t0 = coef( estLpmLin )["age"],
  +     allXVal = xMeanLin["age"], xPos = 1, iPos = 0, model = "lpm" ) )
  + }
         age
  [1,] 42.54
  > if( maxLikLoaded ) {
  +     print( numericGradient( function( x, ... ){ urbinEla( x, ... )$semEla },
  +     t0 = coef( estLpmLin ),
  +     allXVal = xMeanLin, xPos = 3, model = "lpm" ) )
  + }
       (Intercept) kids   age educ
  [1,]           0    0 42.54    0
  > # semi-elasticity of age with standard errors (full covariance matrix)
  > urbinEla( coef( estLpmLin ), xMeanLin, xPos = 3, model = "lpm",
  +   vcov( estLpmLin ) )
   semEla   stdEr 
  -0.3044  0.1114 
  > # semi-elasticity of age with standard errors (only standard errors)
  > urbinEla( coef( estLpmLin )["age"], xMeanLin["age"], xPos = 1, iPos = 0,
  +   model = "lpm", sqrt( diag( vcov( estLpmLin ) ) )["age"] )
   semEla   stdEr 
  -0.3044  0.1114 
  > urbinEla( coef( estLpmLin ), xMeanLin, xPos = 3, model = "lpm",
  +   sqrt( diag( vcov( estLpmLin ) ) ) )
   semEla   stdEr 
  -0.3044  0.1114 
  > 
  > ### quadratic in age
  > estLpmQuad <- lm( lfp ~ kids + age + I(age^2) + educ,
  +   data = Mroz87 )
  > summary( estLpmQuad )
  
  Call:
  lm(formula = lfp ~ kids + age + I(age^2) + educ, data = Mroz87)
  
  Residuals:
     Min     1Q Median     3Q    Max 
  -0.869 -0.531  0.270  0.428  0.742 
  
  Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
  (Intercept) -1.095252   0.518174   -2.11  0.03487 *  
  kids        -0.157282   0.047256   -3.33  0.00092 ***
  age          0.070988   0.024515    2.90  0.00389 ** 
  I(age^2)    -0.000921   0.000287   -3.21  0.00140 ** 
  educ         0.039008   0.007764    5.02  6.3e-07 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 0.482 on 748 degrees of freedom
  Multiple R-squared:  0.0595,	Adjusted R-squared:  0.0545 
  F-statistic: 11.8 on 4 and 748 DF,  p-value: 2.54e-09
  
  > # 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( estLpmQuad )[ c( "age", "I(age^2)" ) ], xMeanQuad[ "age" ],
  +   xPos = c( 1, 2 ), iPos = 0, model = "lpm" )
   semEla   stdEr 
  -0.3122      NA 
  > urbinEla( coef( estLpmQuad ), xMeanQuad, xPos = c( 3, 4 ), model = "lpm" )
   semEla   stdEr 
  -0.3122      NA 
  > # semi-elasticity of age based on numerical derivation
  > 100 * ( predict( estLpmQuad,
  +   newdata = as.data.frame( t( xMeanQuad * c( 1, 1, 1.005, 1.005^2, 1 ) ) ) ) -
  +     predict( estLpmQuad,
  +       newdata = as.data.frame( t( xMeanQuad * c( 1, 1, 0.995, 0.995^2, 1 ) ) ) ) )
        1 
  -0.3122 
  > # partial derivatives of the semi-elasticity wrt the coefficients
  > urbinEla( coef( estLpmQuad )[ c( "age", "I(age^2)" ) ], xMeanQuad[ "age" ],
  +   xPos = c( 1, 2 ), iPos = 0, model = "lpm" )$derivCoef
  [1]   42.54 3618.94
  > urbinEla( coef( estLpmQuad ), xMeanQuad, xPos = c( 3, 4 ), model = "lpm" )$derivCoef
  [1]    0.00    0.00   42.54 3618.94    0.00
  > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEla( x, ... )$semEla },
  +     t0 = coef( estLpmQuad )[ c( "age", "I(age^2)" ) ],
  +     allXVal = xMeanQuad[ "age" ], xPos = c( 1, 2 ), iPos = 0, model = "lpm" ) )
  + }
         age I(age^2)
  [1,] 42.54     3619
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEla( x, ... )$semEla },
  +     t0 = coef( estLpmQuad ),
  +     allXVal = xMeanQuad, xPos = c( 3, 4 ), model = "lpm" ) )
  + }
       (Intercept) kids   age I(age^2) educ
  [1,]           0    0 42.54     3619    0
  > # semi-elasticity of age with standard errors (full covariance matrix)
  > urbinEla( coef( estLpmQuad )[ c( "age", "I(age^2)" ) ], xMeanQuad["age"],
  +   xPos = c( 1, 2 ), iPos = 0, model = "lpm",
  +   vcov( estLpmQuad )[ c( "age", "I(age^2)" ), c( "age", "I(age^2)" ) ] )
   semEla   stdEr 
  -0.3122  0.1107 
  > urbinEla( coef( estLpmQuad ), xMeanQuad, xPos = c( 3, 4 ), model = "lpm",
  +   vcov( estLpmQuad ) )
   semEla   stdEr 
  -0.3122  0.1107 
  > # semi-elasticity of age with standard errors (only standard errors)
  > urbinEla( coef( estLpmQuad )[ c( "age", "I(age^2)" ) ], xMeanQuad[ "age" ],
  +   xPos = c( 1, 2 ), iPos = 0, model = "lpm",
  +   sqrt( diag( vcov( estLpmQuad ) ) )[ c( "age", "I(age^2)" ) ] )
   semEla   stdEr 
  -0.3122  1.4723 
  Warning message:
  In urbinEla(allCoef = coef(estLpmQuad)[c("age", "I(age^2)")], allXVal = xMeanQuad["age"], xPos = c(1, 2), model = "lpm", allCoefVcov = sqrt(diag(vcov(estLpmQuad)))[c("age", "I(age^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 
  > urbinEla( coef( estLpmQuad ), xMeanQuad, xPos = c( 3, 4 ), model = "lpm",
  +   sqrt( diag( vcov( estLpmQuad ) ) ) )
   semEla   stdEr 
  -0.3122  1.4723 
  Warning message:
  In urbinEla(allCoef = coef(estLpmQuad), allXVal = xMeanQuad, xPos = c(3, 4), model = "lpm", allCoefVcov = sqrt(diag(vcov(estLpmQuad)))) :
    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( estLpmQuad )[ c( "age", "I(age^2)" ) ], xMeanQuad["age"],
  +   xPos = c( 1, 2 ), iPos = 0, model = "lpm",
  +   sqrt( diag( vcov( estLpmQuad ) ) )[ c( "age", "I(age^2)" ) ],
  +   xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) )
   semEla   stdEr 
  -0.3122  0.1337 
  > urbinEla( coef( estLpmQuad ), xMeanQuad, xPos = c( 3, 4 ), model = "lpm",
  +   sqrt( diag( vcov( estLpmQuad ) ) ),
  +   xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) )
   semEla   stdEr 
  -0.3122  0.1337 
  > 
  > ### 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
  > estLpmInt <- lm( lfp ~ kids + age30.37 + age38.44 + age53.60 + educ,
  +   data = Mroz87 )
  > summary( estLpmInt )
  
  Call:
  lm(formula = lfp ~ kids + age30.37 + age38.44 + age53.60 + educ, 
      data = Mroz87)
  
  Residuals:
     Min     1Q Median     3Q    Max 
  -0.858 -0.523  0.279  0.437  0.780 
  
  Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
  (Intercept)   0.17366    0.10237    1.70   0.0902 .  
  kids         -0.13594    0.04666   -2.91   0.0037 ** 
  age30.37TRUE  0.05019    0.04773    1.05   0.2934    
  age38.44TRUE  0.04545    0.05130    0.89   0.3759    
  age53.60TRUE -0.19102    0.06105   -3.13   0.0018 ** 
  educ          0.03960    0.00777    5.10  4.4e-07 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 0.483 on 747 degrees of freedom
  Multiple R-squared:  0.0547,	Adjusted R-squared:  0.0484 
  F-statistic: 8.65 on 5 and 747 DF,  p-value: 5.54e-08
  
  > # coefficients of the 'intervals'
  > coefLpmInt <- c( coef( estLpmInt )[3:4], 0, coef( estLpmInt )[5] )
  > # 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] )
  > # mean shares of the 'intervals'
  > xMeanIntShares <- c( xMeanInt[3:4], 1 - sum( xMeanInt[3:5] ), xMeanInt[5] )
  > # semi-elasticity of age without standard errors
  > urbinElaInt( coef( estLpmInt )[3:5], xMeanInt[3:5],
  +   c( 30, 37.5, 44.5, 52.5, 60 ), xPos = c( 1, 2, 0, 3 ), iPos = 0,
  +   model = "lpm" )
   semEla   stdEr 
  -0.3559      NA 
  > urbinElaInt( coef( estLpmInt ), xMeanInt,
  +   c( 30, 37.5, 44.5, 52.5, 60 ), xPos = c( 3, 4, 0, 5 ), model = "lpm" )
   semEla   stdEr 
  -0.3559      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( estLpmInt, newdata = Mroz87Upper ) -
  +     predict( estLpmInt, newdata = Mroz87Lower ) )
  [1] -0.3475
  > Mroz87LowerMean <- Mroz87Lower
  > Mroz87UpperMean <- Mroz87Upper
  > Mroz87LowerMean$kids <- Mroz87UpperMean$kids <- xMeanInt[ "kids" ]
  > Mroz87LowerMean$educ <- Mroz87UpperMean$educ <- xMeanInt[ "educ" ]
  > 10 * mean( predict( estLpmInt, newdata = Mroz87UpperMean ) -
  +     predict( estLpmInt, newdata = Mroz87LowerMean ) )
  [1] -0.3475
  > # partial derivatives of the semi-elasticity wrt the coefficients
  > urbinElaInt( coef( estLpmInt )[3:5], xMeanInt[3:5],
  +   c( 30, 37.5, 44.5, 52.5, 60 ), xPos = c( 1, 2, 0, 3 ), iPos = 0,
  +   model = "lpm" )$derivCoef
  [1] -1.4594 -0.1294  1.4487
  > urbinElaInt( coef( estLpmInt ), xMeanInt,
  +   c( 30, 37.5, 44.5, 52.5, 60 ), xPos = c( 3, 4, 0, 5 ), model = "lpm" )$derivCoef
  [1]  0.0000  0.0000 -1.4594 -0.1294  1.4487  0.0000
  > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinElaInt( x, ... )$semEla },
  +     t0 = coef( estLpmInt )[3:5], allXVal = xMeanInt[3:5],
  +     xBound = c( 30, 37.5, 44.5, 52.5, 60 ), xPos = c( 1, 2, 0, 3 ), iPos = 0,
  +     model = "lpm" ) )
  + }
       age30.37TRUE age38.44TRUE age53.60TRUE
  [1,]       -1.459      -0.1294        1.449
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinElaInt( x, ... )$semEla },
  +     t0 = coef( estLpmInt ), allXVal = xMeanInt,
  +     xBound = c( 30, 37.5, 44.5, 52.5, 60 ), xPos = c( 3, 4, 0, 5 ),
  +     model = "lpm" ) )
  + }
       (Intercept) kids age30.37TRUE age38.44TRUE age53.60TRUE educ
  [1,]           0    0       -1.459      -0.1294        1.449    0
  > # semi-elasticity of age with standard errors (full covariance matrix)
  > vcovLpmInt <- vcov( estLpmInt )
  > vcovLpmInt <- rbind( vcovLpmInt[ 3:4, ], 0, vcovLpmInt[ 5, ] )
  > vcovLpmInt <- cbind( vcovLpmInt[ , 3:4 ], 0, vcovLpmInt[ , 5 ] )
  > urbinElaInt( coef( estLpmInt )[3:5], xMeanInt[3:5],
  +   c( 30, 37.5, 44.5, 52.5, 60 ), xPos = c( 1, 2, 0, 3 ), iPos = 0,
  +   model = "lpm", allCoefVcov = vcov( estLpmInt )[ 3:5, 3:5 ] )
    semEla    stdEr 
  -0.35586  0.09966 
  > urbinElaInt( coef( estLpmInt ), xMeanInt,
  +   c( 30, 37.5, 44.5, 52.5, 60 ), xPos = c( 3, 4, 0, 5 ), model = "lpm",
  +   allCoefVcov = vcov( estLpmInt ) )
    semEla    stdEr 
  -0.35586  0.09966 
  > # semi-elasticity of age with standard errors (only standard errors)
  > urbinElaInt( coef( estLpmInt )[3:5], xMeanInt[3:5],
  +   c( 30, 37.5, 44.5, 52.5, 60 ), xPos = c( 1, 2, 0, 3 ), iPos = 0,
  +   model = "lpm", allCoefVcov = sqrt( diag( vcov( estLpmInt ) ) )[3:5] )
   semEla   stdEr 
  -0.3559  0.1128 
  > urbinElaInt( coef( estLpmInt ), xMeanInt,
  +   c( 30, 37.5, 44.5, 52.5, 60 ), xPos = c( 3, 4, 0, 5 ), model = "lpm",
  +   allCoefVcov = sqrt( diag( vcov( estLpmInt ) ) ) )
   semEla   stdEr 
  -0.3559  0.1128 
  > 
  > 
  > ### 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( estLpmLin )[3], NULL, c( 30, 40 ), c( 50, 60 ), xPos = 1,
  +   iPos = 0, model = "lpm" )
   effect   stdEr 
  -0.1431      NA 
  > urbinEffInt( coef( estLpmLin ), NULL, c( 30, 40 ), c( 50, 60 ), xPos = 3,
  +   model = "lpm" )
   effect   stdEr 
  -0.1431      NA 
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # based on predicted values
  > predict( estLpmLin,
  +   newdata = as.data.frame( t( replace( xMeanLin, 3, 55 ) ) ) ) -
  +   predict( estLpmLin,
  +     newdata = as.data.frame( t( replace( xMeanLin, 3, 35 ) ) ) )
        1 
  -0.1431 
  > # partial derivatives of the semi-elasticity wrt the coefficients
  > urbinEffInt( coef( estLpmLin ), NULL, xPos = 3,
  +   c( 30, 40 ), c( 50, 60 ), model = "lpm" )$derivCoef
  [1]  0  0 20  0
  > # numerically computed partial derivatives of the semi-elasticity wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffInt( x, ... )$effect },
  +     t0 = coef( estLpmLin )[3], allXVal = NULL,
  +     refBound = c( 30, 40 ), intBound = c( 50, 60 ),
  +     xPos = 1, iPos = 0, model = "lpm" ) )
  + }
       age
  [1,]  20
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffInt( x, ... )$effect },
  +     t0 = coef( estLpmLin ), allXVal = NULL,
  +     refBound = c( 30, 40 ), intBound = c( 50, 60 ), xPos = 3, model = "lpm" ) )
  + }
       (Intercept) kids age educ
  [1,]           0    0  20    0
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # (full covariance matrix)
  > urbinEffInt( coef( estLpmLin ), NULL,
  +   c( 30, 40 ), c( 50, 60 ), xPos = 3, model = "lpm",
  +   allCoefVcov = vcov( estLpmLin ) )
    effect    stdEr 
  -0.14310  0.05236 
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # (only standard errors)
  > urbinEffInt( coef( estLpmLin )[3], NULL, c( 30, 40 ), c( 50, 60 ),
  +   xPos = 1, iPos = 0, model = "lpm",
  +   allCoefVcov = sqrt( diag( vcov( estLpmLin ) ) )[3] )
    effect    stdEr 
  -0.14310  0.05236 
  > urbinEffInt( coef( estLpmLin ), NULL, c( 30, 40 ), c( 50, 60 ), xPos = 3,
  +   model = "lpm", allCoefVcov = sqrt( diag( vcov( estLpmLin ) ) ) )
    effect    stdEr 
  -0.14310  0.05236 
  > 
  > 
  > ### 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( estLpmQuad )[3:4], NULL,
  +   c( 30, 40 ), c( 50, 60 ), xPos = 1:2, iPos = 0, model = "lpm" )
   effect   stdEr 
  -0.2375      NA 
  > urbinEffInt( coef( estLpmQuad ), NULL,
  +   c( 30, 40 ), c( 50, 60 ), xPos = 3:4, model = "lpm" )
   effect   stdEr 
  -0.2375      NA 
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # based on predicted values
  > predict( estLpmQuad,
  +   newdata = as.data.frame( t( replace( xMeanQuad, 3:4, c( 55, 55^2 ) ) ) ),
  +   type = "response" ) -
  +   predict( estLpmQuad,
  +     newdata = as.data.frame( t( replace( xMeanQuad, 3:4, c( 35, 35^2 ) ) ) ),
  +     type = "response" )
        1 
  -0.2375 
  > # partial derivatives of the effect wrt the coefficients
  > urbinEffInt( coef( estLpmQuad ), NULL, xPos = c( 3, 4 ),
  +   c( 30, 40 ), c( 50, 60 ), model = "lpm" )$derivCoef
  [1]    0    0   20 1800    0
  > # numerically computed partial derivatives of the effect wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffInt( x, ... )$effect },
  +     t0 = coef( estLpmQuad )[3:4], allXVal = NULL,
  +     refBound = c( 30, 40 ), intBound = c( 50, 60 ), xPos = 1:2, iPos = 0,
  +     model = "lpm" ) )
  + }
       age I(age^2)
  [1,]  20     1800
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffInt( x, ... )$effect },
  +     t0 = coef( estLpmQuad ),
  +     refBound = c( 30, 40 ), intBound = c( 50, 60 ), xPos = 3:4, model = "lpm" ) )
  + }
       (Intercept) kids age I(age^2) educ
  [1,]           0    0  20     1800    0
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # (full covariance matrix)
  > urbinEffInt( coef( estLpmQuad )[3:4], NULL, c( 30, 40 ), c( 50, 60 ),
  +   xPos = 1:2, iPos = 0, model = "lpm",
  +   allCoefVcov = vcov( estLpmQuad )[3:4,3:4] )
    effect    stdEr 
  -0.23747  0.05978 
  > urbinEffInt( coef( estLpmQuad ), NULL, c( 30, 40 ), c( 50, 60 ),
  +   xPos = 3:4, model = "lpm", allCoefVcov = vcov( estLpmQuad ) )
    effect    stdEr 
  -0.23747  0.05978 
  > # effects of age changing from the 30-40 interval to the 50-60 interval
  > # (only standard errors)
  > urbinEffInt( coef( estLpmQuad )[3:4], NULL, c( 30, 40 ), c( 50, 60 ),
  +   xPos = 1:2, iPos = 0, model = "lpm",
  +   allCoefVcov = sqrt( diag( vcov( estLpmQuad ) ) )[3:4] )
   effect   stdEr 
  -0.2375  0.7125 
  Warning message:
  In urbinEffInt(allCoef = coef(estLpmQuad)[3:4], allXVal = NULL, xPos = 1:2, refBound = c(30, 40), intBound = c(50, 60), model = "lpm", allCoefVcov = sqrt(diag(vcov(estLpmQuad)))[3:4], 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 
  > urbinEffInt( coef( estLpmQuad ), NULL, c( 30, 40 ), c( 50, 60 ),
  +   xPos = 3:4, model = "lpm",
  +   allCoefVcov = sqrt( diag( vcov( estLpmQuad ) ) ) )
   effect   stdEr 
  -0.2375  0.7125 
  Warning message:
  In urbinEffInt(allCoef = coef(estLpmQuad), allXVal = NULL, xPos = 3:4, refBound = c(30, 40), intBound = c(50, 60), model = "lpm", allCoefVcov = sqrt(diag(vcov(estLpmQuad)))) :
    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( estLpmQuad ), NULL, xPos = c( 3, 4 ),
  +   c( 30, 40 ), c( 50, 60 ), model = "lpm",
  +   allCoefVcov = sqrt( diag( vcov( estLpmQuad ) ) ),
  +   xMeanSd = c( mean( Mroz87$age ), sd( Mroz87$age ) ) )
    effect    stdEr 
  -0.23747  0.06992 
  > 
  > ### 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( estLpmInt )[3:5], xMeanInt[3:5], 1:3, iPos = 0,
  +   c( -1, -1, 1, 0 ), model = "lpm" )
   effect   stdEr 
  -0.2393      NA 
  > urbinEffCat( coef( estLpmInt ), xMeanInt, 3:5, c( -1, -1, 1, 0 ),
  +   model = "lpm" )
   effect   stdEr 
  -0.2393      NA 
  > # effects calculated based on predicted values
  > names( xMeanInt ) <- sub( "TRUE", "", names( coef( estLpmInt ) ) )
  > 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( estLpmInt, newdata = df53.60 ) -
  +   sum( Mroz87$age30.37 ) / sum( Mroz87$age30.37 + Mroz87$age38.44 ) *
  +   predict( estLpmInt, newdata = df30.37 ) -
  +   sum( Mroz87$age38.44 ) / sum( Mroz87$age30.37 + Mroz87$age38.44 ) *
  +   predict( estLpmInt, newdata = df38.44 )
        1 
  -0.2393 
  > # partial derivatives of the effect wrt the coefficients
  > urbinEffCat( coef( estLpmInt )[3:5], xMeanInt[3:5], 1:3, iPos = 0,
  +   c( -1, -1, 1, 0 ), model = "lpm" )$derivCoef
  [1] -0.5878 -0.4122  1.0000
  > urbinEffCat( coef( estLpmInt ), xMeanInt, 3:5, c( -1, -1, 1, 0 ),
  +   model = "lpm" )$derivCoef
  [1]  0.0000  0.0000 -0.5878 -0.4122  1.0000  0.0000
  > # numerically computed partial derivatives of the effect wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffCat( x, ... )$effect },
  +     t0 = coef( estLpmInt )[3:5], xPos = 1:3, iPos = 0,
  +     allXVal = xMeanInt[3:5], xGroups = c( -1, -1, 1, 0 ), model = "lpm" ) )
  + }
       age30.37TRUE age38.44TRUE age53.60TRUE
  [1,]      -0.5878      -0.4122            1
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffCat( x, ... )$effect },
  +     t0 = coef( estLpmInt ), xPos = 3:5,
  +     allXVal = xMeanInt, xGroups = c( -1, -1, 1, 0 ), model = "lpm" ) )
  + }
       (Intercept) kids age30.37TRUE age38.44TRUE age53.60TRUE educ
  [1,]           0    0      -0.5878      -0.4122            1    0
  > # with full covariance matrix
  > urbinEffCat( coef( estLpmInt )[3:5], xMeanInt[3:5], 1:3, iPos = 0,
  +   c( -1, -1, 1, 0 ), model = "lpm",
  +   allCoefVcov = vcov( estLpmInt )[3:5, 3:5] )
    effect    stdEr 
  -0.23926  0.06451 
  > urbinEffCat( coef( estLpmInt ), xMeanInt, 3:5, c( -1, -1, 1, 0 ),
  +   model = "lpm", allCoefVcov = vcov( estLpmInt ) )
    effect    stdEr 
  -0.23926  0.06451 
  > # with standard errors only
  > urbinEffCat( coef( estLpmInt )[3:5], xMeanInt[3:5], 1:3, iPos = 0,
  +   c( -1, -1, 1, 0 ), model = "lpm",
  +   allCoefVcov = sqrt( diag( vcov( estLpmInt ) ) )[3:5] )
    effect    stdEr 
  -0.23926  0.07044 
  > urbinEffCat( coef( estLpmInt ), xMeanInt, 3:5, c( -1, -1, 1, 0 ),
  +   model = "lpm", allCoefVcov = sqrt( diag( vcov( estLpmInt ) ) ) )
    effect    stdEr 
  -0.23926  0.07044 
  > ### effects of age changing from the 53-60 category to the 38-52 category
  > # without standard errors
  > urbinEffCat( coef( estLpmInt )[3:5], xMeanInt[3:5], 1:3, iPos = 0,
  +   c( 0, 1, -1, 1 ), model = "lpm" )
  effect  stdEr 
  0.2108     NA 
  > urbinEffCat( coef( estLpmInt ), xMeanInt, 3:5, c( 0, 1, -1, 1 ),
  +   model = "lpm" )
  effect  stdEr 
  0.2108     NA 
  > # effects calculated based on predicted values
  > sum( Mroz87$age38.44 ) / sum( Mroz87$age38.44 + Mroz87$age45.52 ) *
  +   predict( estLpmInt, newdata = df38.44 ) +
  +   sum( Mroz87$age45.52 ) / sum( Mroz87$age38.44 + Mroz87$age45.52 ) *
  +   predict( estLpmInt, newdata = df45.52 ) -
  +   predict( estLpmInt, newdata = df53.60 )
       1 
  0.2108 
  > # partial derivatives of the effect wrt the coefficients
  > urbinEffCat( coef( estLpmInt )[3:5], xMeanInt[3:5], 1:3, iPos = 0,
  +   c( 0, 1, -1, 1 ), model = "lpm" )$derivCoef
  [1]  0.0000  0.4346 -1.0000
  > urbinEffCat( coef( estLpmInt ), xMeanInt, 3:5, c( 0, 1, -1, 1 ),
  +   model = "lpm" )$derivCoef
  [1]  0.0000  0.0000  0.0000  0.4346 -1.0000  0.0000
  > # numerically computed partial derivatives of the effect wrt the coefficients
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffCat( x, ... )$effect },
  +     t0 = coef( estLpmInt )[3:5], xPos = 1:3, iPos = 0,
  +     allXVal = xMeanInt[3:5], xGroups = c( 0, 1, -1, 1 ), model = "lpm" ) )
  + }
       age30.37TRUE age38.44TRUE age53.60TRUE
  [1,]            0       0.4346           -1
  > if( maxLikLoaded ) {
  +   print( numericGradient( function( x, ... ){ urbinEffCat( x, ... )$effect },
  +     t0 = coef( estLpmInt ), xPos = 3:5,
  +     allXVal = xMeanInt, xGroups = c( 0, 1, -1, 1 ), model = "lpm" ) )
  + }
       (Intercept) kids age30.37TRUE age38.44TRUE age53.60TRUE educ
  [1,]           0    0            0       0.4346           -1    0
  > # with full covariance matrix
  > urbinEffCat( coef( estLpmInt )[3:5], xMeanInt[3:5], 1:3, iPos = 0,
  +   c( 0, 1, -1, 1 ), model = "lpm",
  +   allCoefVcov = vcov( estLpmInt )[3:5,3:5] )
  effect  stdEr 
  0.2108 0.0599 
  > urbinEffCat( coef( estLpmInt ), xMeanInt, 3:5, c( 0, 1, -1, 1 ),
  +   model = "lpm", allCoefVcov = vcov( estLpmInt ) )
  effect  stdEr 
  0.2108 0.0599 
  > # with standard errors only
  > urbinEffCat( coef( estLpmInt )[3:5], xMeanInt[3:5], 1:3, iPos = 0,
  +   c( 0, 1, -1, 1 ), model = "lpm",
  +   allCoefVcov = sqrt( diag( vcov( estLpmInt ) ) )[3:5] )
   effect   stdEr 
  0.21077 0.06499 
  > urbinEffCat( coef( estLpmInt ), xMeanInt, 3:5, c( 0, 1, -1, 1 ),
  +   model = "lpm", allCoefVcov = sqrt( diag( vcov( estLpmInt ) ) ) )
   effect   stdEr 
  0.21077 0.06499 
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.01    0.20    2.18 
Running the tests in 'tests/probit.R' failed.
Complete output:
  > 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 
     2.68    0.26    2.89