Date: 2012-07-05 15:17 Sender: Tomas Krehlikok, I added truncated and censored means and put it as a function:
pred <- function(model, data, type = "fitted") {
coefficients <- summary(model)$estimate[1:(length(summary(model)$estimate[,1])-1),1]
names <- names(coefficients)
indices <- unlist(lapply(names[-1], function(x) which(colnames(data)==x)))
if (length(indices)==(length(coefficients)-1)) {
data_pred <- cbind(rep(1, length(data[,1])), do.call(cbind, lapply(indices, function(x) data[,x])))
fitted <- (data_pred %*% coefficients)[,1]
} else {
warning("The dataset does not contain some of the variables.")
}
sigma <- summary(model.tobit)$estimate[length(summary(model.tobit)$estimate[,1]),1]
if (type == "fitted") {
fitted
} else if ((type == "censored") || (type == "truncated")) {
if (model$left == -Inf) {
print("Right censored/truncated mean:")
right <- model$right
truncated <- fitted - sigma*(lambda((fitted - right)/sigma))
prob_trunc <- 1 - pnorm((fitted - right)/sigma)
censored <- right*dnorm((fitted - right)/sigma) + prob_trunc*truncated
} else if (model$right == Inf) {
print("Left censored/truncated mean:")
left <- model$left
truncated <- fitted + sigma*(lambda((fitted - left)/sigma))
prob_trunc <- pnorm((fitted - left)/sigma)
censored <- left*dnorm((fitted - left)/sigma) + prob_trunc*truncated
} else {
print("Interval censored/truncated mean:")
left <- model$left
right <- model$right
truncated <- fitted - sigma*((dnorm((left-fitted)/sigma)-dnorm((right-fitted)/sigma))/(pnorm((left-fitted)/sigma)-pnorm((right-fitted)/sigma)))
prob_trunc <- pnorm((fitted - right)/sigma) - pnorm((fitted - left)/sigma)
censored <- left*dnorm((fitted - left)/sigma) + prob_trunc*truncated + right*dnorm((fitted - right)/sigma)
}
if (type == "censored") {
censored
} else {
truncated
}
}
}
lambda <- function(x) {
lambda <- dnorm(x)/pnorm(x)
lambda
} |