R Under development (unstable) (2024-09-30 r87204 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library( "urbin" ) > maxLikLoaded <- require( "maxLik" ) Loading required package: maxLik Loading required package: miscTools Please cite the 'maxLik' package as: Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1. If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site: https://r-forge.r-project.org/projects/maxlik/ > 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 1.64 0.14 1.76