# load freight data
data(freight)
# Compute Standard Poisson estimates
glm_model <- glm(broken ~ transfers, data=freight,
family=poisson, na.action=na.exclude) # beta estimates
print("The standard Poisson estimates for the beta vector are")
print(coef(glm_model))
# Compute COM-Poisson estimates (under constant dispersion model)
start.time <- Sys.time()
cmp_model = cmp(formula = broken ~ transfers, data=freight)
print("The COM-Poisson estimates for the beta vector are")
print(coef(cmp_model))
print("The COM-Poisson estimate for the dispersion parameter nu is")
print(nu(cmp_model))
# Compute associated standard errors for constant COM-Poisson estimates
print("The associated standard errors for the betas in the constant dispersion case are")
print(sdev(cmp_model))
# Perform likelihood ratio test for dispersion parameter
# Test for dispersion equal or not equal to 1 (ie performing Poisson vs COM-Poisson regression)
freight.chisq <- chisq(cmp_model)
freight.pval <- pval(cmp_model)
print(sprintf("The likelihood ratio chi-squared test statistic is %0.5f
and associated p-value (testing Poisson vs CMP regression) is %0.5f",
freight.chisq, freight.pval))
# Compute constant COM-Poisson leverage
freight.lev <- leverage(cmp_model)
print("The leverage of the points is")
print(freight.lev)
# Compute constant COM-Poisson deviances
# commented-out for speed
# freight.CMPDev <- deviance(cmp_model)
# print("The approximate constant dispersion standardized CMP Deviance is")
# print(freight.CMPDev)
# Compute fitted values
freight.fitted = predict(cmp_model, newdata=freight)
print("The CMP fitted values are")
print(freight.fitted)
# Compute residual values
freight.constantCMPresids <- residuals(cmp_model)
print("The CMP residuals are")
print(freight.constantCMPresids)
# Compute MSE
freight.constantCMP.MSE <- mean(freight.constantCMPresids^2)
print("The MSE for the constant CMP regression is")
print(freight.constantCMP.MSE)
# Compute predictions on new data
new_data = data.frame(transfers=(0:10))
freight.predicted = predict(cmp_model, newdata=new_data)
plot(0:10, freight.predicted, type="l",
xlab="number of transfers", ylab="predicted number broken")
# Compute parametric bootstrap results and use them to generate
# 0.95 confidence intervals for parameters
# commented-out for speed
# freight.CMPParamBoot <- parametric_bootstrap(cmp_model, n=1000)
# print(apply(freight.CMPParamBoot,2,quantile,c(0.025,0.975)))
Run the code above in your browser using DataLab