# Poisson-Tweedie: properties and regression models for count data ----- # Author: Wagner Hugo Bonat LEG/IMADA ---------------------------------- # Script 5 - Customer profile ------------------------------------------ # Date: 13/07/2016 ----------------------------------------------------- rm(list=ls()) # Loading extra packages require(mcglm) require(xtable) # Loading data set data <- read.table("PTWDataset4.txt") # Linear predictor form <- ncust ~ nhu + aid + aha + dnc + ds # Poisson model fit0 <- glm(form, family = poisson, data = data) # Initial values ------------------------------------------------------- ini <- list() ini$regression <- list("R1" = coef(fit0)) ini$tau <- list("R1" = c(1)) ini$rho = c(0) # Neyman-Type A ini$power <- list("R1" = c(1)) fit1 <- mcglm(c(form), list(mc_id(data)), link = "log", power_fixed = TRUE, variance = "poisson_tweedie", control_ini = ini, data = data) summary(fit1) # Negative binomial ini$power <- list("R1" = c(2)) fit2 <- mcglm(c(form), list(mc_id(data)), link = "log", power_fixed = TRUE, variance = "poisson_tweedie", control_ini = ini, control_algorithm = list(tunning = 0.75), data = data) summary(fit2) # Poisson-inverse Gaussian ini$power <- list("R1" = c(3)) fit3 <- mcglm(c(form), list(mc_id(data)), link = "log", power_fixed = TRUE, variance = "poisson_tweedie", control_ini = ini, control_algorithm = list(tunning = 0.75), data = data) summary(fit3) # Measures of goodness-of-fit rbind(gof(fit1),gof(fit2),gof(fit3)) # Regression coefficients xtable(round(cbind(coef(fit0),coef(fit1, type = "beta")$Estimates, coef(fit2, type = "beta")$Estimates, coef(fit3, type = "beta")$Estimates),3), digits = 3) # Standard errors round(cbind(sqrt(diag(vcov(fit0))), coef(fit1, type = "beta",std.error = TRUE)$Std.error, coef(fit2, type = "beta",std.error = TRUE)$Std.error, coef(fit3, type = "beta",std.error = TRUE)$Std.error),3) round(rbind(coef(fit1, type = "tau",std.error = TRUE)[,1:2], coef(fit2, type = "tau",std.error = TRUE)[,1:2], coef(fit3, type = "tau",std.error = TRUE)[,1:2]),3) # END ------------------------------------------------------------------ save.image("Script5.RData")