2

To analyse right-censored data, I'm using the tobit regression package called censReg in R. When assuming a fixed effect model, this works fine.

library(censReg)
library(tidyverse)
x=rep(c(0,2,4,6,8),20)
beta0=rep(rnorm(20,0,1),each=5)
ind=rep(1:20,each=5)
beta=2
y=x*beta + beta0 +rnorm(100,0,0.5)
data=data.frame(x=x,y=y,ind=ind) %>% mutate(y_cens=ifelse(y>12,12,y))

fit = censReg(y_cens ~ x,left=-Inf,right=12, data = data)

fit

Call: censReg(formula = y_cens ~ x, left = -Inf, right = 12, data = data)

Coefficients: (Intercept) x logSigma -0.35715
2.00235 0.01105

But when assuming a random intercept (as I would do in other R packages), the function does not seem to work anymore.

fit = censReg(y_cens ~ x + (1|ind),left=-Inf,right=12, data = data)

fit

Call: censReg(formula = y_cens ~ x + (1 | ind), left = -Inf, right = 12, data = data)

Coefficients: NULL

However,in the vignette (see here), the authors mention some methods to solve random-effect model, but they do not explain how to write it.

Anthony
  • 441
  • 2
  • 10

1 Answers1

2

You cannot specify a random intercept with the (1|ind) call in censReg() . You need to create a new data object with pdata.frame() from the plm package, giving the column names that contain identifiers of participants (ind) and the time variable (x). This new data object will be a list, not a data frame. Becausethe first two columns in the data object created with pdata.frame are factors, you need to create a new column with numeric data (x1 below). Then you specify the method as method = "BHHH". This is all explained in the vigniette :-)

library(censReg)
library(plm)
pData <- pdata.frame(data, c( "ind", "x" ) )
pData$x1 <- as.numeric(pData$x)
fit <- censReg(y_cens ~ x1, data = pData, method = "BHHH", left=-Inf, right=12)
summary(fit)
Call:
censReg(formula = y_cens ~ x1, left = -Inf, right = 12, data = pData, 
    method = "BHHH")

Observations:
         Total  Left-censored     Uncensored Right-censored 
           100              0             68             32


Coefficients:
            Estimate Std. error t value Pr(> t)

(Intercept)  -3.8876     0.3581 -10.855  <2e-16 ***
x1            4.0168     0.1353  29.684  <2e-16 ***
logSigma      0.2311     0.1196   1.932  0.0534 .




Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


BHHH maximisation, 32 iterations
Return code 8: successive function values within relative tolerance limit (reltol)
Log-likelihood: -119.9258 on 3 Df
Uki Buki
  • 101