# Poisson-Tweedie: properties and regression models for count data ----- # Author: Wagner Hugo Bonat LEG/IMADA ---------------------------------- # Script 2 - Morbidity by respiratory disease in Curitiba/PR - Brazil--- # Date: 11/07/2016 ----------------------------------------------------- rm(list=ls()) # Loading extra packages require(mcglm) require(MASS) require(Matrix) # Loading data set data <- read.table("PTWDataset1.txt", header= TRUE) # Linear predictor form <- MORB04 ~ cos(2*pi*MONTH/12) + sin(2*pi*MONTH/12) + MAXIMA + MINIMA + PRECIPITACAO # Fitting Poisson regression model fit0 <- glm(form, family = poisson, data = data, offset = log(data$POP04)) # Fitting Poisson-Tweedie regression model fit1 <- mcglm(c(form), list(mc_id(data)), link = "log", variance = "poisson_tweedie", power_fixed = FALSE, data = data, offset = list(c(log(data$POP04))), control_algorithm = list(verbose = FALSE, max_iter = 100, tunning = 0.15, correct = TRUE)) summary(fit1) # Computing conditional standard error for dispersion parameter std.phi <- mc_conditional_test(fit1, parameters = c("power11","tau11"), test = 2, fixed = 1) # Table 2 -------------------------------------------------------------- cbind(coef(fit0), coef(fit1, type = "beta")$Estimates, coef(fit1, type = "beta")$Estimates/coef(fit0)) cbind(sqrt(diag(vcov(fit0))), coef(fit1, type = "beta", std.error = TRUE)$Std.error, coef(fit1, type = "beta", std.error = TRUE)$Std.error /sqrt(diag(vcov(fit0)))) rbind(coef(fit1, type = "power", std.error = TRUE)[,1:2], std.phi[1:2]) # Figure 6 ------------------------------------------------------------- # Fitted values and standard errors X <- model.matrix(form, data = data) eta <- X%*%coef(fit1, type = "beta")$Estimates std.error <- sqrt(diag(as.matrix(X%*%vcov(fit1)[1:6,1:6]%*%t(X)))) pred <- exp(eta + log(data$POP04)) lim.min <- exp(eta - qnorm(0.975)*std.error + log(data$POP04)) lim.max <- exp(eta + qnorm(0.975)*std.error + log(data$POP04)) # Layout of the Figure mat <- matrix(NA,ncol=3,nrow=2) mat[1,1:3] <- 1 mat[2,1] <- 2 mat[2,2] <- 3 mat[2,3] <- 4 layout(mat) # Plot par(mar=c(2.8, 2.8, 2, 0.5), mgp = c(1.6, 0.6, 0)) with(data, plot(MORB04 ~ TREND, type = "p", ylab = "Morbidity", xlab = "Month/Year", xaxt = "n", main = "A", col = "gray", pch = 20)) labels <- data$DATA posi <- seq(from = 1, to=length(labels), by=5) axis(1, at = data$TREND[posi], labels = labels[posi], las = 1, cex.axis=0.8) lines(data$TREND, pred, lty = 2) lines(data$TREND, lim.min) lines(data$TREND, lim.max) with(data, plot(MORB04 ~ PRECIPITACAO, ylab = "Morbidity", pch = 20, xlab = "Precipitation", main = "B", col = "gray")) abline(coef(lm(MORB04 ~ PRECIPITACAO, data = data))) with(data, plot(MORB04 ~ MAXIMA, ylab = "Morbidity", pch = 20, xlab = "Maxima", main = "C", col = "gray")) abline(coef(lm(MORB04 ~ MAXIMA, data = data))) with(data, plot(MORB04 ~ MINIMA, ylab = "Morbidity", pch = 20, xlab = "Minima", main = "D", col = "gray")) abline(coef(lm(MORB04 ~ MINIMA, data = data))) # END ------------------------------------------------------------------ save.image("Script2.RData")