Weighted censored regression

Question

As continuation of this question

I have a model with a continuous predictor and factors. I want the fit of the the model (especially the beta to the continuous predictor) to fit best when 2 factors are "on". I tried to use weights. The results in the toy-example are alright but I can not improve my real world problem.

The example is:

library(AER)
N = 10
f = rep(c("s1","s2","s3","s4","s5","s6","s7","s8"),N)
fcoeff = rep(c(-1,-2,-3,-4,-3,-5,-10,-5),N)
set.seed(100) 
x = rnorm(8*N)+1
beta = 5
epsilon = rnorm(8*N,sd = sqrt(1/5))
y.star = x*beta+fcoeff+epsilon ## latent response
y = y.star 
y[y<0] <- 0 ## censored response

my.weights = rep(c(1,1,1,2,2,1,1,1),N)

fit <- tobit(y~0+x+f, weights = my.weights)
summary(fit)

coef(fit) 

Giving more weight (2 instead of 1) at $s4$ and $s5$ I try to force the model to fit these data points more exactly as errors here are worse than errors elsewhere.

Is this a valid approach? Can this be done differently?

EDIT after Achim Zeileis' answer: I tried to add the interaction with $s4$ and $s5$ using 2 sets of factors:

f = rep(c("s1","s2","s3","s4","s5","s6","s7","s8"),N)
f2 = rep(c("s1","s2","s3","s4","s5","s6","s7","s8"),N)
f[f %in% c("s4","s5")] <- "no.inter"
f2[f2 %in% c("s1","s2","s3","s6","s7","s8")] <- "rest"

fit <- tobit(y~0+x+x*f2+f)

The results is

coef(fit)

         x     f2rest       f2s4       f2s5        fs1        fs2        fs3        fs6        fs7        fs8     x:f2s4     x:f2s5 
 4.9662343 -4.4759248 -3.7568885 -4.4251267  3.4118712  2.4866271  1.2519721 -0.5528713 -5.4638955  0.0000000 -0.2309364  1.5865312 

Thus I get the correct overall beta of 5. A level shift of approx -4.5 and shifts of $s_i,i \in \{1,2,3,6,7,8\}$ that look good with the last one being $0$ thus approximatively giving the shift of $-5$ in connection with the factor f2rest. Finally there are the interactions with s4 and s5. This looks like the solution. I try to model interaction at those special factors only. I have not seen anyone doing this before. Is it correct that I need to define those 2 sets of factors?

EDIT: I have mixed something up. The above model is singular. I am confused. I don't know how I made the first fit work.

Answer 1

The approach is certainly ad hoc so that valid might be too strong. Whether or not it yields useful results, depends on what you really want to do with that model.

One aspect that you will lose are the standard errors. tobit() and survreg() treat the weights as case weights (like most other maximum likelihood regression functions in R). This means that doubling the weights doubles the information and decreases the variances by 1/2. For example:

fit1 <- tobit(y ~ 0 + x + f, weights = rep(1, 8 * N))
fit10 <- tobit(y ~ 0 + x + f, weights = rep(10, 8 * N))
all.equal(coef(fit10), coef(fit1))
## [1] TRUE
all.equal(vcov(fit10), vcov(fit1))
## [1] "Mean relative difference: 9"

Thus, the coefficients are not affected but the estimated variances of the coefficient estimates are deacreased by a factor of 1/10.

Note that this is different from most implementations of least-squares linear regression where weights are typically just used as proportionality weights but not as case weights. This is also the case in R, see e.g., this discussion: https://mailman.stat.ethz.ch/pipermail/r-help/2012-February/302563.html

An alternative approach to adapting the slope for the continuous regressor x would be to include interactions for all or for selected levels of the factor f. And then you can test whether slopes are significantly different etc.

EDIT: Follow-up on the idea with the interaction. If you want different intercepts with respect to all 8 groups and different slopes with respect to only s4/s5 vs. the rest, you can define to factors f and f2 almost as you did above.

f <- factor(f)
f2 <- factor(f %in% c("s4", "s5"))

The code is somewhat simpler and it is cleaner to use factor variables rather than character variables as you did. Then the desired interaction model is:

fit <- tobit(y ~ 0 + f + x + f2:x)
coef(fit)
##        fs1        fs2        fs3        fs4        fs5 
## -1.0677546 -1.9939367 -3.2296259 -3.9268500 -2.7746633 
##        fs6        fs7        fs8          x   x:f2TRUE 
## -5.0351431 -9.9448424 -4.4811849  4.9686998 -0.1258431 

This shows first the eight different intercepts, then the slope of x in the group f2 = FALSE, and finally the difference in slopes wrt x in the group f2 = TRUE (i.e., s4/s5). So in the latter group, the slope is 4.84 rather than 4.97. This interpretation of the model formula is maybe easiest to see when you look at the model matrix:

head(model.matrix(fit))
##   fs1 fs2 fs3 fs4 fs5 fs6 fs7 fs8         x x:f2TRUE
## 1   1   0   0   0   0   0   0   0 0.4978076 0.000000
## 2   0   1   0   0   0   0   0   0 1.1315312 0.000000
## 3   0   0   1   0   0   0   0   0 0.9210829 0.000000
## 4   0   0   0   1   0   0   0   0 1.8867848 1.886785
## 5   0   0   0   0   1   0   0   0 1.1169713 1.116971
## 6   0   0   0   0   0   1   0   0 1.3186301 0.000000

The first nine columns should be clear. And the tenth column is the same as x if f2 = TRUE (i.e., f in s4/s5) and zero otherwise.

For more details on R's formula language and the interpretation of interactions, you can look at most textbooks that discuss linear models with R or tutorials on the web etc.

Answer 2

I have another viewpoint, based on the preference of weighting the regression when two factors are 'turned on'. In essence, you have a preferential weighting of the regression line in the presence of selected factors. Now, I quoted a passage from an introduction to Bayesian regression to quote:

When we want show the linear fit from a Bayesian model, instead of showing only estimate, we can draw a range of lines, with each one representing a different estimate of the model parameters. As the number of datapoints increases, the lines begin to overlap because there is less uncertainty in the model parameters.

Sounds familiar? Also, with respect to the interpretation of parameters to quote:

The result of performing Bayesian Linear Regression is a distribution of possible model parameters based on the data and the prior. This allows us to quantify our uncertainty about the model: if we have fewer data points, the posterior distribution will be more spread out. As the amount of data points increases, the likelihood washes out the prior, and in the case of infinite data, the outputs for the parameters converge to the values obtained from OLS.

Bottom line, my suggestion adopt a more formalistic Bayesian linear regression approach.