将 vcov 矩阵添加到 AER::tobit

如何解决将 vcov 矩阵添加到 AER::tobit

我想为 vcov 提供一个 AER::tobit 矩阵。

如果我查看 https://rdrr.io/cran/AER/src/R/tobit.R，我会看到这个选项：

vcov.tobit <- function(object,...) {
  vc <- NextMethod()
  if(is.null(colnames(vc))) {
    nam <- names(object$coefficients)    
    nam <- if(length(nam) == ncol(vc)) {
      nam
    } else if(length(nam) == ncol(vc) - 1L) {
      c(nam,"Log(scale)")
    } else {
      c(nam,names(object$scale))
    }
    colnames(vc) <- rownames(vc) <- nam
  }
  return(vc)  
}

但是，如果执行以下操作：

library(AER)
s1.lm <- lm(taxrate ~ votewon + industry + size + urbanisation + vote,data=DF)
yhat <- fitted(s1.lm)
s2.tobit <- AER::tobit(sales ~ yhat + industry + size + urbanisation + vote,left=12,right=33,data=DF)
DF <- cbind(DF,yhat)
vcov_mat = vcov2sls(s1.lm,s2.tobit,DF)

现在我想添加 vcov

s2.tobit <- AER::tobit(sales ~ yhat + industry + size + urbanisation + vote,vcov.tobit=vcov_mat,data=DF)

或

s2.tobit_corrected <- AER::vcov.tobit(sales ~ yhat + industry + size + urbanisation + vote,vcov=vcov_mat,data=DF)

我明白了：

Error in survreg.control(...) : unused argument (vcov.tobit = vcov_mat)

还有：

Error: 'vcov.tobit' is not an exported object from 'namespace:AER'

我应该如何添加这个 vcov matix？

数据

vcov2sls <- function(s1,s2,data,type=2) {
  ## get y names
  y1.nm <- gsub(".*=\\s(.*)(?=\\s~).*","\\1",deparse(s1$call)[1],perl=TRUE)
  y2.nm <- "sales"
  ## auxilliary model matrix
  X <- cbind(`(Intercept)`=1,data[,y1.nm,F],model.matrix(s2)[,-(1:2)])
  ## get y
  y <- data[,y2.nm] 
  ## betas second stage
  b <- s2$coefficients
  ## calculate corrected sums of squares
  sse <- sum((y - b %*% t(X))^2)
  rmse <- sqrt(mean(s2$residuals^2))  ## RMSE 2nd stage
  V0 <- vcov(s2)  ## biased vcov 2nd stage
  dof <- s2$df.residual  ## degrees of freedom 2nd stage
  ## calculate corrected RMSE
  rmse.c <- sqrt(sse/dof)
  ## calculate corrected vcov
  V <- (rmse.c/rmse)^2 * V0
  return(V)
}

DF <- structure(list(country = c("C","C","J","B","F","E","D","I","H","G","A","F"),year = c(2005,2010,2005,2005),sales = c(15.48,12.39,3.72,23.61,4,31.87,25.33,7.64,-0.26,2.9,15.48,2.9),industry = c("D",urbanisation = c("B","B"),size = c(1,1,5,2,3,5),base_rate = c(14L,14L,19L,30L,20L,29L,24L,17L,33L,23L,20L),taxrate = c(12L,12L,21L,18L,32L,22L,25L,vote = c(0,1),votewon = c(0,1)),class = "data.frame",row.names = c(NA,-100L))

## convert variables to factors beforehand
DF[c(1,6,9,10)] <- lapply(DF[c(1,10)],factor)