SCM

[#2126] predict for censReg

Date:
2012-07-05 12:38
Priority:
3
State:
Open
Submitted by:
Tomas Krehlik (tomaskrehlik)
Assigned to:
Arne Henningsen (arne)
Product:
censReg
Operating System:
All
Status:
Summary:
predict for censReg

Detailed description
Hello,

I pretty miss some predict functionality in censReg, that would at provide fitted values. And to be constructive, I am appending code, which seems to be pretty robust for usage :) (do not know how slow it is, but I use it regularly and serves well)

(just replace model.tobit by generic model and data3 by name of the provided dataset and probably some functionality which would care about if there is a constatnt in the model - honestly, I am not sure if censored models without constant make any sense )

# Get coefficients
coefficients <- summary(model.tobit)$estimate[1:(length(summary(model.tobit)$estimate[,1])-1),1]
# Get names of the coefficients
names <- names(coefficients)
# Find indices of the names
indices <- unlist(lapply(names[-1], function(x) which(colnames(data3)==x)))
# Check if they exist in the dataset (if they were found)
if (length(indices)==(length(coefficients)-1)) {
# Add constant term and construct dataset
data <- cbind(rep(1, length(data3[,1])), do.call(cbind, lapply(indices, function(x) data3[,x])))
# Make the predicted values
fitted <- (data %*% coefficients)[,1]
} else {
warning("The dataset does not contain some of the variables.")
}

Best,
TK.

Comments:

Message  ↓
Date: 2012-07-06 03:52
Sender: Arne Henningsen

Dear Tomas
Thanks for submitting the Feature Request and for submitting code for this. I will try to implement it ASAP but I am probably not able to do this before my summer holidays.
/Arne

Date: 2012-07-05 15:17
Sender: Tomas Krehlik

ok, 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
}

Attached Files:

Changes

Field Old Value Date By
assigned_tonone2012-07-06 03:49arne
ProductNone2012-07-06 03:49arne
Operating SystemNone2012-07-06 03:49arne
Thanks to:
Vienna University of Economics and Business Powered By FusionForge