R version 4.4.0 RC (2024-04-16 r86451 ucrt) -- "Puppy Cup" 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( "censReg" ) 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/ Please cite the 'censReg' package as: Henningsen, Arne (2017). censReg: Censored Regression (Tobit) Models. R package version 0.5. http://CRAN.R-Project.org/package=censReg. If you have questions, suggestions, or comments regarding the 'censReg' package, please use a forum or 'tracker' at the R-Forge site of the 'sampleSelection' project: https://r-forge.r-project.org/projects/sampleselection/ > library( "plm" ) > > # load outputs that were previously produced by this script > saved <- new.env() > load( "censRegPanelTest.RData.save", envir = saved ) > > options( digits = 5 ) > > printAll <- function( objName, what = "print" ) { + cat( "Comparing new object '", objName, "' to previously saved object...", + sep = "" ) + x <- get( objName ) + if( !exists( objName, envir = saved, inherits = FALSE ) ) { + cat( " previously saved object not found\n" ) + } else { + xSaved <- get( objName, envir = saved, inherits = FALSE ) + if( !isTRUE( all.equal( class( x ), class( xSaved ) ) ) ) { + cat( " different classes:\n" ) + cat( "new:\n" ) + print( class( x ) ) + cat( "saved:\n" ) + print( class( xSaved ) ) + } else if( !isTRUE( all.equal( names( x ), names( xSaved ) ) ) ) { + cat( " different names:\n" ) + cat( "new:\n" ) + print( names( x ) ) + cat( "saved:\n" ) + print( names( xSaved ) ) + } else { + cat( "\n" ) + } + for( n in names( x ) ) { + if( ! n %in% c( "code", "gradient", "iterations", "last.step", + "message" ) ) { + cat( " comparing component '", n, "' ...", sep = "" ) + if( n == "vcov" ) { + tol <- 5e-1 + } else if( n == "estimate" ) { + tol <- 5e-2 + } else { + tol <- 5e-3 + } + testRes <- all.equal( x[[ n ]], xSaved[[ n ]], tol = tol ) + if( isTRUE( testRes ) ) { + cat( " OK\n" ) + } else { + cat( " different\n" ) + print( testRes ) + cat( "new:\n" ) + print( x[[ n ]] ) + cat( "saved:\n" ) + try( print( xSaved[[ n ]] ) ) + } + } + } + } + + for( mName in c( "Coef", "CoefNoLs", "Vcov", "VcovNoLs", + "CoefSum", "CoefSumNoLs", "LogLik", "Nobs", "ExtractAIC" ) ) { + if( mName == "Vcov" & objName == "randEffOnlyInt" ) { + next + } + cat( " comparing method '", mName, "' ...", sep = "" ) + tol <- 5e-3 + if( mName == "Coef" ) { + xm <- coef( x ) + tol <- 5e-2 + } else if( mName == "CoefNoLs" ) { + xm <- coef( x, logSigma = FALSE ) + tol <- 5e-2 + } else if( mName == "Vcov" ) { + xm <- vcov( x ) + tol <- 5e-1 + } else if( mName == "VcovNoLs" ) { + xm <- vcov( x, logSigma = FALSE ) + tol <- 5e-1 + } else if( mName == "CoefSum" ) { + xm <- coef( summary( x ) ) + tol <- 5e-2 + } else if( mName == "CoefSumNoLs" ) { + xm <- coef( summary( x ), logSigma = FALSE ) + tol <- 5e-2 + } else if( mName == "LogLik" ) { + xm <- logLik( x ) + } else if( mName == "Nobs" ) { + xm <- nobs( x ) + } else if( mName == "ExtractAIC" ) { + xm <- extractAIC( x ) + } else { + stop( "unknown value of 'mName': ", mName ) + } + methodObjName <- paste0( objName, mName ) + if( !exists( methodObjName, envir = saved, inherits = FALSE ) ) { + cat( " previously saved object not found\n" ) + } else { + xmSaved <- get( methodObjName, envir = saved, inherits = FALSE ) + testRes <- all.equal( xm, xmSaved, tol = tol ) + if( isTRUE( testRes ) ) { + cat( " OK\n" ) + } else { + cat( " different\n" ) + print( testRes ) + cat( "new:\n" ) + print( xm ) + cat( "saved:\n" ) + print( xmSaved ) + } + } + # assign to parent frame so that it will be included in the saved workspace + assign( methodObjName, xm, envir = parent.frame() ) + } + + if( what %in% c( "print", "methods", "all" ) ) { + print( x, digits = 1 ) + print( x, logSigma = FALSE , digits = 1 ) + print( maxLik:::summary.maxLik( x ), digits = 1 ) + print( summary( x ), digits = 1 ) + print( summary( x ), logSigma = FALSE , digits = 1 ) + } + if( what %in% c( "methods", "all" ) ) { + print( round( coef( x ), 2 ) ) + print( round( coef( x, logSigma = FALSE ), 2 ) ) + print( round( vcov( x ), 2 ) ) + print( round( vcov( x, logSigma = FALSE ), 2 ) ) + print( round( coef( summary( x ) ), 2 ) ) + print( round( coef( summary( x ), logSigma = FALSE ), 2 ) ) + try( margEff( x ) ) + print( logLik( x ) ) + print( nobs( x ) ) + print( extractAIC( x ) ) + } + + if( what == "all" ) { + for( n in names( x ) ) { + cat( "$", n, "\n", sep = "" ) + if( n %in% c( "estimate", "gradientObs" ) ) { + print( round( x[[ n ]], 2 ) ) + } else if( n %in% c( "hessian" ) ) { + print( round( x[[ n ]], 1 ) ) + } else if( n %in% c( "gradient" ) ) { + } else if( ! n %in% c( "last.step" ) ) { + print( x[[ n ]] ) + } + cat( "\n" ) + } + cat( "class\n" ) + print( class( x ) ) + } + } > > nId <- 15 > nTime <- 4 > > set.seed( 123 ) > pData <- data.frame( + id = rep( paste( "F", 1:nId, sep = "_" ), each = nTime ), + time = rep( 1980 + 1:nTime, nId ) ) > pData$ui <- rep( rnorm( nId ), each = nTime ) > pData$x1 <- rnorm( nId * nTime ) > pData$x2 <- runif( nId * nTime ) > pData$ys <- -1 + pData$ui + 2 * pData$x1 + 3 * pData$x2 + rnorm( nId * nTime ) > pData$y <- ifelse( pData$ys > 0, pData$ys, 0 ) > nData <- pData # save data set without information on panel structure > pData <- pdata.frame( pData, c( "id", "time" ), stringsAsFactors = FALSE ) > > > ## Newton-Raphson method > randEff <- censReg( y ~ x1 + x2, data = pData ) > printAll( "randEff" ) Comparing new object 'randEff' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'Vcov' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK Call: censReg(formula = y ~ x1 + x2, data = pData) Coefficients: (Intercept) x1 x2 logSigmaMu logSigmaNu -0.37 1.68 2.24 -0.13 -0.01 Call: censReg(formula = y ~ x1 + x2, data = pData) Coefficients: (Intercept) x1 x2 sigmaMu sigmaNu -0.4 1.7 2.2 0.9 1.0 -------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-Likelihood: -73.199 5 free parameters Estimates: Estimate Std. error t value Pr(> t) (Intercept) -0.37 0.47 -0.8 0.4 x1 1.68 0.21 8.0 1e-15 *** x2 2.24 0.67 3.3 9e-04 *** logSigmaMu -0.13 0.26 -0.5 0.6 logSigmaNu -0.01 0.13 -0.1 0.9 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -------------------------------------------- Call: censReg(formula = y ~ x1 + x2, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.37 0.47 -0.8 0.4 x1 1.68 0.21 8.0 1e-15 *** x2 2.24 0.67 3.3 9e-04 *** logSigmaMu -0.13 0.26 -0.5 0.6 logSigmaNu -0.01 0.13 -0.1 0.9 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -73.199 on 5 Df Call: censReg(formula = y ~ x1 + x2, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.4 0.5 -0.8 0.4 x1 1.7 0.2 8.0 1e-15 *** x2 2.2 0.7 3.3 9e-04 *** sigmaMu 0.9 0.2 3.9 1e-04 *** sigmaNu 1.0 0.1 7.7 1e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -73.199 on 5 Df > try( margEff( randEff ) ) Error in margEff.censReg(randEff) : the margEff() method for objects of class 'censReg' can not yet be used for panel data models > # only intercept > randEffOnlyInt <- censReg( y ~ 1, data = pData ) > printAll( "randEffOnlyInt", what = "diff" ) Comparing new object 'randEffOnlyInt' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK Warning messages: 1: In sqrt(diag(vc)) : NaNs produced 2: In sqrt(diag(vc)) : NaNs produced 3: In sqrt(diag(vc)) : NaNs produced 4: In sqrt(diag(vc)) : NaNs produced > # no intercept > randEffNoInt <- censReg( y ~ x1 -1, data = pData ) > printAll( "randEffNoInt" ) Comparing new object 'randEffNoInt' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'Vcov' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK Call: censReg(formula = y ~ x1 - 1, data = pData) Coefficients: x1 logSigmaMu logSigmaNu 1.9 0.1 0.2 Call: censReg(formula = y ~ x1 - 1, data = pData) Coefficients: x1 sigmaMu sigmaNu 2 1 1 -------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-Likelihood: -81.186 3 free parameters Estimates: Estimate Std. error t value Pr(> t) x1 1.9 0.3 7.3 2e-13 *** logSigmaMu 0.1 0.2 0.5 0.6 logSigmaNu 0.2 0.1 1.6 0.1 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -------------------------------------------- Call: censReg(formula = y ~ x1 - 1, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) x1 1.9 0.3 7.3 2e-13 *** logSigmaMu 0.1 0.2 0.5 0.6 logSigmaNu 0.2 0.1 1.6 0.1 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -81.186 on 3 Df Call: censReg(formula = y ~ x1 - 1, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) x1 1.9 0.3 7 2e-13 *** sigmaMu 1.1 0.3 4 6e-05 *** sigmaNu 1.2 0.2 8 2e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -81.186 on 3 Df > # neither intercept nor explanatory variables > try( censReg( y ~ -1, data = pData ) ) Error in censReg(y ~ -1, data = pData) : argument 'formula' seems to have neither an intercept nor any explanatory variables > > ## BHHH method > randEffBhhh <- censReg( y ~ x1 + x2, data = pData, method = "BHHH" ) > printAll( "randEffBhhh" ) Comparing new object 'randEffBhhh' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'Vcov' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK Call: censReg(formula = y ~ x1 + x2, data = pData, method = "BHHH") Coefficients: (Intercept) x1 x2 logSigmaMu logSigmaNu -0.37 1.68 2.24 -0.13 -0.01 Call: censReg(formula = y ~ x1 + x2, data = pData, method = "BHHH") Coefficients: (Intercept) x1 x2 sigmaMu sigmaNu -0.4 1.7 2.2 0.9 1.0 -------------------------------------------- Maximum Likelihood estimation BHHH maximisation, 14 iterations Return code 8: successive function values within relative tolerance limit (reltol) Log-Likelihood: -73.199 5 free parameters Estimates: Estimate Std. error t value Pr(> t) (Intercept) -0.37 0.56 -0.7 0.510 x1 1.68 0.29 5.7 1e-08 *** x2 2.24 0.73 3.1 0.002 ** logSigmaMu -0.13 0.29 -0.4 0.660 logSigmaNu -0.01 0.14 -0.1 0.930 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -------------------------------------------- Call: censReg(formula = y ~ x1 + x2, data = pData, method = "BHHH") Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.37 0.56 -0.7 0.510 x1 1.68 0.29 5.7 1e-08 *** x2 2.24 0.73 3.1 0.002 ** logSigmaMu -0.13 0.29 -0.4 0.660 logSigmaNu -0.01 0.14 -0.1 0.930 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 BHHH maximisation, 14 iterations Return code 8: successive function values within relative tolerance limit (reltol) Log-likelihood: -73.199 on 5 Df Call: censReg(formula = y ~ x1 + x2, data = pData, method = "BHHH") Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.4 0.6 -0.7 0.510 x1 1.7 0.3 5.7 1e-08 *** x2 2.2 0.7 3.1 0.002 ** sigmaMu 0.9 0.3 3.4 7e-04 *** sigmaNu 1.0 0.1 7.1 1e-12 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 BHHH maximisation, 14 iterations Return code 8: successive function values within relative tolerance limit (reltol) Log-likelihood: -73.199 on 5 Df > > > ## BFGS method (optim) > randEffBfgs <- censReg( y ~ x1 + x2, data = pData, method = "BFGS" ) > printAll( "randEffBfgs" ) Comparing new object 'randEffBfgs' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'constraints' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'Vcov' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK Call: censReg(formula = y ~ x1 + x2, data = pData, method = "BFGS") Coefficients: (Intercept) x1 x2 logSigmaMu logSigmaNu -0.37 1.68 2.24 -0.13 -0.01 Call: censReg(formula = y ~ x1 + x2, data = pData, method = "BFGS") Coefficients: (Intercept) x1 x2 sigmaMu sigmaNu -0.4 1.7 2.2 0.9 1.0 -------------------------------------------- Maximum Likelihood estimation BFGS maximization, 25 iterations Return code 0: successful convergence Log-Likelihood: -73.199 5 free parameters Estimates: Estimate Std. error t value Pr(> t) (Intercept) -0.37 0.47 -0.8 0.4 x1 1.68 0.21 8.0 1e-15 *** x2 2.24 0.67 3.3 9e-04 *** logSigmaMu -0.13 0.26 -0.5 0.6 logSigmaNu -0.01 0.13 -0.1 0.9 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -------------------------------------------- Call: censReg(formula = y ~ x1 + x2, data = pData, method = "BFGS") Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.37 0.47 -0.8 0.4 x1 1.68 0.21 8.0 1e-15 *** x2 2.24 0.67 3.3 9e-04 *** logSigmaMu -0.13 0.26 -0.5 0.6 logSigmaNu -0.01 0.13 -0.1 0.9 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 BFGS maximization, 25 iterations Return code 0: successful convergence Log-likelihood: -73.199 on 5 Df Call: censReg(formula = y ~ x1 + x2, data = pData, method = "BFGS") Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.4 0.5 -0.8 0.4 x1 1.7 0.2 8.0 1e-15 *** x2 2.2 0.7 3.3 9e-04 *** sigmaMu 0.9 0.2 3.9 1e-04 *** sigmaNu 1.0 0.1 7.7 1e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 BFGS maximization, 25 iterations Return code 0: successful convergence Log-likelihood: -73.199 on 5 Df > > > ## BFGS method (R) > randEffBfgsr <- censReg( y ~ x1 + x2, data = pData, method = "BFGSR" ) > printAll( "randEffBfgsr", what = "none" ) Comparing new object 'randEffBfgsr' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'Vcov' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK > > > ## BHHH with starting values > randEffBhhhStart <- censReg( y ~ x1 + x2, data = pData, method = "BHHH", + start = c( -0.4, 1.7, 2.2, -0.1, -0.01 ) ) > printAll( "randEffBhhhStart" ) Comparing new object 'randEffBhhhStart' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'Vcov' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK Call: censReg(formula = y ~ x1 + x2, data = pData, start = c(-0.4, 1.7, 2.2, -0.1, -0.01), method = "BHHH") Coefficients: (Intercept) x1 x2 logSigmaMu logSigmaNu -0.37 1.68 2.24 -0.13 -0.01 Call: censReg(formula = y ~ x1 + x2, data = pData, start = c(-0.4, 1.7, 2.2, -0.1, -0.01), method = "BHHH") Coefficients: (Intercept) x1 x2 sigmaMu sigmaNu -0.4 1.7 2.2 0.9 1.0 -------------------------------------------- Maximum Likelihood estimation BHHH maximisation, 7 iterations Return code 8: successive function values within relative tolerance limit (reltol) Log-Likelihood: -73.199 5 free parameters Estimates: Estimate Std. error t value Pr(> t) (Intercept) -0.37 0.56 -0.7 0.510 x1 1.68 0.29 5.7 1e-08 *** x2 2.24 0.73 3.1 0.002 ** logSigmaMu -0.13 0.30 -0.4 0.661 logSigmaNu -0.01 0.14 -0.1 0.929 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -------------------------------------------- Call: censReg(formula = y ~ x1 + x2, data = pData, start = c(-0.4, 1.7, 2.2, -0.1, -0.01), method = "BHHH") Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.37 0.56 -0.7 0.510 x1 1.68 0.29 5.7 1e-08 *** x2 2.24 0.73 3.1 0.002 ** logSigmaMu -0.13 0.30 -0.4 0.661 logSigmaNu -0.01 0.14 -0.1 0.929 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 BHHH maximisation, 7 iterations Return code 8: successive function values within relative tolerance limit (reltol) Log-likelihood: -73.199 on 5 Df Call: censReg(formula = y ~ x1 + x2, data = pData, start = c(-0.4, 1.7, 2.2, -0.1, -0.01), method = "BHHH") Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.4 0.6 -0.7 0.510 x1 1.7 0.3 5.7 1e-08 *** x2 2.2 0.7 3.1 0.002 ** sigmaMu 0.9 0.3 3.4 7e-04 *** sigmaNu 1.0 0.1 7.1 1e-12 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 BHHH maximisation, 7 iterations Return code 8: successive function values within relative tolerance limit (reltol) Log-likelihood: -73.199 on 5 Df > > > ## left-censoring at 5 > pData$yAdd <- pData$y + 5 > randEffAdd <- censReg( yAdd ~ x1 + x2, data = pData, left = 5 ) > printAll( "randEffAdd" ) Comparing new object 'randEffAdd' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'Vcov' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK Call: censReg(formula = yAdd ~ x1 + x2, left = 5, data = pData) Coefficients: (Intercept) x1 x2 logSigmaMu logSigmaNu 4.63 1.68 2.24 -0.13 -0.01 Call: censReg(formula = yAdd ~ x1 + x2, left = 5, data = pData) Coefficients: (Intercept) x1 x2 sigmaMu sigmaNu 4.6 1.7 2.2 0.9 1.0 -------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-Likelihood: -73.199 5 free parameters Estimates: Estimate Std. error t value Pr(> t) (Intercept) 4.63 0.47 9.8 <2e-16 *** x1 1.68 0.21 8.0 1e-15 *** x2 2.24 0.67 3.3 9e-04 *** logSigmaMu -0.13 0.26 -0.5 0.6 logSigmaNu -0.01 0.13 -0.1 0.9 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -------------------------------------------- Call: censReg(formula = yAdd ~ x1 + x2, left = 5, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) 4.63 0.47 9.8 <2e-16 *** x1 1.68 0.21 8.0 1e-15 *** x2 2.24 0.67 3.3 9e-04 *** logSigmaMu -0.13 0.26 -0.5 0.6 logSigmaNu -0.01 0.13 -0.1 0.9 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -73.199 on 5 Df Call: censReg(formula = yAdd ~ x1 + x2, left = 5, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 20 40 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) 4.6 0.5 10 <2e-16 *** x1 1.7 0.2 8 1e-15 *** x2 2.2 0.7 3 9e-04 *** sigmaMu 0.9 0.2 4 1e-04 *** sigmaNu 1.0 0.1 8 1e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -73.199 on 5 Df > > > ## right-censoring > pData$yNeg <- - pData$y > randEffNeg <- censReg( yNeg ~ x1 + x2, data = pData, + left = -Inf, right = 0 ) > printAll( "randEffNeg" ) Comparing new object 'randEffNeg' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'Vcov' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK Call: censReg(formula = yNeg ~ x1 + x2, left = -Inf, right = 0, data = pData) Coefficients: (Intercept) x1 x2 logSigmaMu logSigmaNu 0.37 -1.68 -2.24 -0.13 -0.01 Call: censReg(formula = yNeg ~ x1 + x2, left = -Inf, right = 0, data = pData) Coefficients: (Intercept) x1 x2 sigmaMu sigmaNu 0.4 -1.7 -2.2 0.9 1.0 -------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-Likelihood: -73.199 5 free parameters Estimates: Estimate Std. error t value Pr(> t) (Intercept) 0.37 0.47 0.8 0.4 x1 -1.68 0.21 -8.0 1e-15 *** x2 -2.24 0.67 -3.3 9e-04 *** logSigmaMu -0.13 0.26 -0.5 0.6 logSigmaNu -0.01 0.13 -0.1 0.9 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -------------------------------------------- Call: censReg(formula = yNeg ~ x1 + x2, left = -Inf, right = 0, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 0 40 20 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) 0.37 0.47 0.8 0.4 x1 -1.68 0.21 -8.0 1e-15 *** x2 -2.24 0.67 -3.3 9e-04 *** logSigmaMu -0.13 0.26 -0.5 0.6 logSigmaNu -0.01 0.13 -0.1 0.9 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -73.199 on 5 Df Call: censReg(formula = yNeg ~ x1 + x2, left = -Inf, right = 0, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 0 40 20 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) 0.4 0.5 0.8 0.4 x1 -1.7 0.2 -8.0 1e-15 *** x2 -2.2 0.7 -3.3 9e-04 *** sigmaMu 0.9 0.2 3.9 1e-04 *** sigmaNu 1.0 0.1 7.7 1e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -73.199 on 5 Df > > > ## right-censoring at -5 > pData$yAddNeg <- - pData$yAdd > randEffAddNeg <- censReg( yAddNeg ~ x1 + x2, data = pData, + left = -Inf, right = -5 ) > printAll( "randEffAddNeg" ) Comparing new object 'randEffAddNeg' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'Vcov' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK Call: censReg(formula = yAddNeg ~ x1 + x2, left = -Inf, right = -5, data = pData) Coefficients: (Intercept) x1 x2 logSigmaMu logSigmaNu -4.63 -1.68 -2.24 -0.13 -0.01 Call: censReg(formula = yAddNeg ~ x1 + x2, left = -Inf, right = -5, data = pData) Coefficients: (Intercept) x1 x2 sigmaMu sigmaNu -4.6 -1.7 -2.2 0.9 1.0 -------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-Likelihood: -73.199 5 free parameters Estimates: Estimate Std. error t value Pr(> t) (Intercept) -4.63 0.47 -9.8 <2e-16 *** x1 -1.68 0.21 -8.0 1e-15 *** x2 -2.24 0.67 -3.3 9e-04 *** logSigmaMu -0.13 0.26 -0.5 0.6 logSigmaNu -0.01 0.13 -0.1 0.9 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -------------------------------------------- Call: censReg(formula = yAddNeg ~ x1 + x2, left = -Inf, right = -5, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 0 40 20 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -4.63 0.47 -9.8 <2e-16 *** x1 -1.68 0.21 -8.0 1e-15 *** x2 -2.24 0.67 -3.3 9e-04 *** logSigmaMu -0.13 0.26 -0.5 0.6 logSigmaNu -0.01 0.13 -0.1 0.9 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -73.199 on 5 Df Call: censReg(formula = yAddNeg ~ x1 + x2, left = -Inf, right = -5, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 0 40 20 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -4.6 0.5 -10 <2e-16 *** x1 -1.7 0.2 -8 1e-15 *** x2 -2.2 0.7 -3 9e-04 *** sigmaMu 0.9 0.2 4 1e-04 *** sigmaNu 1.0 0.1 8 1e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -73.199 on 5 Df > > > ## both right and left censoring > pData$yBoth <- ifelse( pData$y < 3, pData$y, 3 ) > randEffBoth <- censReg( yBoth ~ x1 + x2, data = pData, + left = 0, right = 3 ) > printAll( "randEffBoth" ) Comparing new object 'randEffBoth' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'Vcov' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK Call: censReg(formula = yBoth ~ x1 + x2, left = 0, right = 3, data = pData) Coefficients: (Intercept) x1 x2 logSigmaMu logSigmaNu -0.236 1.893 1.972 0.002 0.053 Call: censReg(formula = yBoth ~ x1 + x2, left = 0, right = 3, data = pData) Coefficients: (Intercept) x1 x2 sigmaMu sigmaNu -0.2 1.9 2.0 1.0 1.1 -------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 6 iterations Return code 1: gradient close to zero (gradtol) Log-Likelihood: -64.313 5 free parameters Estimates: Estimate Std. error t value Pr(> t) (Intercept) -0.236 0.549 -0.4 0.67 x1 1.893 0.301 6.3 3e-10 *** x2 1.972 0.819 2.4 0.02 * logSigmaMu 0.002 0.278 0.0 0.99 logSigmaNu 0.053 0.163 0.3 0.75 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -------------------------------------------- Call: censReg(formula = yBoth ~ x1 + x2, left = 0, right = 3, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 20 28 12 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.236 0.549 -0.4 0.67 x1 1.893 0.301 6.3 3e-10 *** x2 1.972 0.819 2.4 0.02 * logSigmaMu 0.002 0.278 0.0 0.99 logSigmaNu 0.053 0.163 0.3 0.75 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 6 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -64.313 on 5 Df Call: censReg(formula = yBoth ~ x1 + x2, left = 0, right = 3, data = pData) Observations: Total Left-censored Uncensored Right-censored 60 20 28 12 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.2 0.5 -0.4 0.67 x1 1.9 0.3 6.3 3e-10 *** x2 2.0 0.8 2.4 0.02 * sigmaMu 1.0 0.3 3.6 3e-04 *** sigmaNu 1.1 0.2 6.1 9e-10 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 6 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -64.313 on 5 Df > > > ## re-order observations/individuals > set.seed( 234 ) > perm <- sample( nId ) > nData2 <- nData > nData2$id <- NA > for( i in 1:nId ) { + nData2$id[ nData$id == paste( "F", i, sep = "_" ) ] <- + paste( "G", perm[ i ], sep = "_" ) + } > pData2 <- pdata.frame( nData2, c( "id", "time" ), stringsAsFactors = FALSE ) > randEff2 <- censReg( y ~ x1 + x2, data = pData2 ) > all.equal( randEff2[ -c(3,5,6,7,9,11,15) ], + randEff[ -c(3,5,6,7,9,11,15) ], tolerance = 1e-2 ) [1] TRUE > > # check if the order of observations/individuals influences the likelihood values > d1c1 <- censReg( y ~ x1 + x2, data = pData, start = coef(randEff), + iterlim = 0 ) > all.equal( d1c1[-c(5,6,7,9,12,15,19)], randEff[-c(5,6,7,9,12,15,19)] ) [1] TRUE > round( d1c1$maximum - randEff$maximum, 12 ) [1] 0 > > d2c2 <- censReg( y ~ x1 + x2, data = pData2, start = coef(randEff2), + iterlim = 0 ) > all.equal( d2c2[-c(5,6,7,9,12,15,19)], randEff2[-c(5,6,7,9,12,15,19)] ) [1] TRUE > round( d2c2$maximum - randEff2$maximum, 12 ) [1] 0 > > d1c2 <- censReg( y ~ x1 + x2, data = pData, + start = coef(randEff2), iterlim = 0 ) > round( d2c2$maximum - d1c2$maximum, 12 ) [1] 0 > round( d2c2$gradient - d1c2$gradient, 12 ) (Intercept) x1 x2 logSigmaMu logSigmaNu 0 0 0 0 0 > > d2c1 <- censReg( y ~ x1 + x2, data = pData2, + start = coef(randEff), iterlim = 0 ) > round( d1c1$maximum - d2c1$maximum, 12 ) [1] 0 > round( d1c1$gradient - d2c1$gradient, 12 ) (Intercept) x1 x2 logSigmaMu logSigmaNu 0 0 0 0 0 > > round( d2c2$maximum - d2c1$maximum, 3 ) [1] 0 > round( d1c1$maximum - d1c2$maximum, 3 ) [1] 0 > > d1cS <- censReg( y ~ x1 + x2, data = pData, + start = randEff$start, iterlim = 0 ) > d2cS <- censReg( y ~ x1 + x2, data = pData2, + start = randEff$start, iterlim = 0 ) > round( d1cS$maximum - d2cS$maximum, 12 ) [1] 0 > round( d1cS$gradient - d2cS$gradient, 12 ) (Intercept) x1 x2 logSigmaMu logSigmaNu 0 0 0 0 0 > > > ## unbalanced panel data > nDataUnb <- nData[ -c( 2, 5, 6, 8 ), ] > pDataUnb <- pdata.frame( nDataUnb, c( "id", "time" ), stringsAsFactors = FALSE ) > randEffUnb <- censReg( y ~ x1 + x2, data = pDataUnb ) > printAll( "randEffUnb" ) Comparing new object 'randEffUnb' to previously saved object... comparing component 'maximum' ... OK comparing component 'estimate' ... OK comparing component 'hessian' ... OK comparing component 'fixed' ... OK comparing component 'type' ... OK comparing component 'gradientObs' ... OK comparing component 'control' ... OK comparing component 'objectiveFn' ... OK comparing component 'xMean' ... OK comparing component 'call' ... OK comparing component 'terms' ... OK comparing component 'nObs' ... OK comparing component 'df.residual' ... OK comparing component 'start' ... OK comparing component 'left' ... OK comparing component 'right' ... OK comparing method 'Coef' ... OK comparing method 'CoefNoLs' ... OK comparing method 'Vcov' ... OK comparing method 'VcovNoLs' ... OK comparing method 'CoefSum' ... OK comparing method 'CoefSumNoLs' ... OK comparing method 'LogLik' ... OK comparing method 'Nobs' ... OK comparing method 'ExtractAIC' ... OK Call: censReg(formula = y ~ x1 + x2, data = pDataUnb) Coefficients: (Intercept) x1 x2 logSigmaMu logSigmaNu -0.225 1.641 2.112 -0.167 -0.001 Call: censReg(formula = y ~ x1 + x2, data = pDataUnb) Coefficients: (Intercept) x1 x2 sigmaMu sigmaNu -0.2 1.6 2.1 0.8 1.0 -------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-Likelihood: -71.193 5 free parameters Estimates: Estimate Std. error t value Pr(> t) (Intercept) -0.225 0.473 -0.5 0.634 x1 1.641 0.211 7.8 8e-15 *** x2 2.112 0.685 3.1 0.002 ** logSigmaMu -0.167 0.271 -0.6 0.538 logSigmaNu -0.001 0.132 0.0 0.993 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -------------------------------------------- Call: censReg(formula = y ~ x1 + x2, data = pDataUnb) Observations: Total Left-censored Uncensored Right-censored 56 17 39 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.225 0.473 -0.5 0.634 x1 1.641 0.211 7.8 8e-15 *** x2 2.112 0.685 3.1 0.002 ** logSigmaMu -0.167 0.271 -0.6 0.538 logSigmaNu -0.001 0.132 0.0 0.993 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -71.193 on 5 Df Call: censReg(formula = y ~ x1 + x2, data = pDataUnb) Observations: Total Left-censored Uncensored Right-censored 56 17 39 0 Coefficients: Estimate Std. error t value Pr(> t) (Intercept) -0.2 0.5 -0.5 0.634 x1 1.6 0.2 7.8 8e-15 *** x2 2.1 0.7 3.1 0.002 ** sigmaMu 0.8 0.2 3.7 2e-04 *** sigmaNu 1.0 0.1 7.6 4e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero (gradtol) Log-likelihood: -71.193 on 5 Df > > > ## NAs in data > pDataNa <- pData > obsNa <- which( ! rownames( pData ) %in% rownames( pDataUnb ) ) > pDataNa$y[ obsNa[ 1:2 ] ] <- NA > pDataNa$x1[ obsNa[ 3 ] ] <- NA > pDataNa$x2[ obsNa[ c( 1, 2, 4 ) ] ] <- NA > randEffNa <- censReg( y ~ x1 + x2, data = pDataNa ) > all.equal( randEffNa[ -15 ], randEffUnb[ -15 ] ) [1] TRUE > > > # returning log-likelihood contributions only (no estimations) > logLikRandEff <- censReg( y ~ x1 + x2, data = pData, start = coef( randEff ), + logLikOnly = TRUE ) > print( logLikRandEff, digits = 1 ) [1] -4 -4 -7 -5 -3 -5 -6 -2 -4 -4 -7 -6 -4 -6 -6 attr(,"gradient") [,1] [,2] [,3] [,4] [,5] [1,] -0.85 0.23 -0.07 -0.1 -1.0 [2,] -1.70 0.05 -1.19 1.7 -1.0 [3,] 1.75 0.33 1.73 1.8 1.8 [4,] 0.16 -0.20 -0.43 -0.7 -1.4 [5,] 0.13 0.98 0.69 -0.7 -0.6 [6,] 0.33 -0.31 0.38 -0.5 -1.6 [7,] -0.08 -3.19 -0.24 -0.7 6.0 [8,] -0.27 0.38 0.06 -0.5 -1.1 [9,] 1.15 0.59 0.30 0.1 -1.6 [10,] -0.44 -0.41 -0.15 -0.5 -1.9 [11,] -0.03 1.45 -0.64 -0.9 2.7 [12,] 1.21 0.54 0.61 0.2 -1.7 [13,] 0.54 -0.23 0.33 -0.4 -1.4 [14,] -1.81 -1.31 -0.92 2.0 2.6 [15,] -0.08 1.09 -0.46 -0.8 0.2 > all.equal( sum( logLikRandEff ), c( logLik( randEff ) ) ) [1] TRUE > logLikStart <- censReg( y ~ x1 + x2, data = pData, + start = c( -0.4, 1.7, 2.2, -0.1, -0.01 ), logLikOnly = TRUE ) > print( round( c( logLikStart ), 3 ) ) [1] -4.224 -3.591 -7.110 -5.316 -3.389 -5.360 -6.349 -1.662 -4.434 -3.774 [11] -6.788 -5.699 -4.385 -5.593 -5.548 > print( round( attr( logLikStart, "gradient" ), 2 ) ) [,1] [,2] [,3] [,4] [,5] [1,] -0.79 0.28 -0.04 -0.18 -1.03 [2,] -1.58 0.09 -1.10 1.49 -1.06 [3,] 1.72 0.26 1.72 1.91 1.84 [4,] 0.18 -0.21 -0.41 -0.61 -1.53 [5,] 0.16 0.97 0.71 -0.66 -0.51 [6,] 0.36 -0.39 0.42 -0.43 -1.58 [7,] -0.05 -3.20 -0.23 -0.70 6.13 [8,] -0.22 0.36 0.08 -0.57 -1.07 [9,] 1.13 0.55 0.30 0.15 -1.61 [10,] -0.35 -0.39 -0.11 -0.52 -1.91 [11,] -0.01 1.42 -0.62 -0.97 2.68 [12,] 1.16 0.50 0.60 0.20 -1.75 [13,] 0.60 -0.27 0.37 -0.30 -1.38 [14,] -1.71 -1.20 -0.85 1.82 2.50 [15,] -0.03 1.06 -0.42 -0.92 0.20 > > > # save all objectives that were produced in this script > # (in order to compare them with objects created by this script in the future) > rm( saved ) > save.image( "censRegPanelTest.RData" ) > > > proc.time() user system elapsed 6.17 0.18 6.37