# Example: Left Truncated outcome 
library(micsr)
data("charitable")

summary(charitable)

charitable$logdon <- log(charitable$donation) - log(25)

char_form <- logdon ~ log(donparents) + log(income) +
  education + religion + married + south
ch_ml <- tobit1(char_form, data = charitable)

summary(ch_ml)

ch_twostep <- tobit1(char_form, data = charitable, method = "twostep")
ch_scls <- tobit1(char_form, data = charitable, method = "trimmed")
ch_ols <- tobit1(char_form, data = charitable, method = "lm")


library(tidyverse)
library(broom)
models <- list(
  OLS      = ch_ols,
  TwoStep  = ch_twostep,
  ML  = ch_ml,
  SCLS  = ch_scls
)

tidy_models <- map_dfr(
  names(models),
  ~ tidy(models[[.x]]) %>% mutate(model = .x)
)

tidy_models <- tidy_models %>%
  filter(term != "(Intercept)", !is.na(estimate), !is.na(std.error))

ggplot(tidy_models, aes(x = term, y = estimate, color = model)) +
  geom_point(
    position = position_dodge(width = 0.6),
    size = 3
  ) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  coord_flip() +
  theme_minimal() +
  labs(
    title = "Coefficient Estimates Across Models",
    x = "",
    y = "Estimated Coefficient",
    color = "Model"
  )

cmtest(ch_ml, test = "normality")
cmtest(ch_ml, test = "heterosc") 

cmtest(ch_ml, test = "skewness")
cmtest(ch_ml, test = "kurtosis")




# Right-truncated outcome
food <- as.data.frame(micsr.data::food)
food <- food[food$year == 1980,]
head(food)
help(food)
summary(food)

plot(food$income, food$food)
abline(h=13030)
sum(food$food==17670)
mean(food$food==17670)

food_tobit <- tobit1(log(food) ~ log(income) + log(hsize) + midage,
                     data = food, #subset = year == 1980,
                     left = -Inf, right = log(13030))
food_tobit 
summary(food_tobit)
##########
# Two-sided tobit model 
portfolio <- as.data.frame(micsr.data::portfolio)
names(portfolio)
summary(portfolio$share)

portfolio$agesq  <- portfolio$age^2/100
prec_ml <- tobit1(share ~ uncert + expinc + networth + 
                    age + agesq + female,
                  left = 0, right = 1, data = portfolio)

prec_ml

prec_ht <- crch::crch(share ~ uncert + expinc + networth +
                        age + agesq + female | networth +
                        age + agesq + female, left = 0, right = 1,
                      data = portfolio)

summary(prec_ht)$coefficients

summary(prec_ht)$loglik
summary(prec_ml)$logLik


prec_ht2 <- crch::crch(share ~ uncert + expinc + networth +
                        age + agesq + female | networth, left = 0, right = 1,
                      data = portfolio)


AIC(prec_ml)
AIC(prec_ht2)
AIC(prec_ht)