Untitled page
PoissonRegression.R
############### Poisson regression
require(ggplot2) # For plotting
require(msm)
library(MASS) # [Link]
library(pscl) # zeroinfl + vuong
library(dplyr)
df<-[Link]("/Users/apple/Downloads/[Link]",header=TRUE)
#Divide the data set into two: 75% train, 25% test
smp_size <- floor(0.75 * nrow(df))
## set the seed to make your partition reproducible
[Link](123)
[Link] <- sample(1:nrow(df), size = smp_size)
train_df <- df[[Link], ]
test_df <- df[-[Link], ]
with(train_df, tapply(travels, lo_card, function(x) {
sprintf("M (SD) = %1.2f (%1.2f)", mean(x), sd(x))
}))
ggplot(train_df, aes(travels, fill = lo_card)) +
geom_histogram(binwidth=.5, position="dodge")
# --- Fit models on TRAIN ----
m_pois <- glm(travels ~ disc_ons + lo_card, family = poisson(link = "log"), data = train_df)
summary(m_pois)
[Link]<-step(m_pois, direction = "backward")
summary([Link])
##################### Issues
100*sum(train_df$travels == 0)/nrow(train_df)
# Negative Binomial Regression
m_nb <- [Link](travels ~ disc_ons + lo_card, data = train_df)
# ZIP: count ~ X | zero-inflation ~ Z
m_zip <- zeroinfl(travels ~ disc_ons + lo_card | disc_ons + lo_card,
dist = "poisson", link = "logit", data = train_df)
#summary(m_pois)
summary(m_nb)
summary(m_zip)
# --- Predict on TEST ----
pred_pois <- predict(m_pois, newdata = test_df, type = "response")
pred_nb <- predict(m_nb, newdata = test_df, type = "response")
pred_zip <- predict(m_zip, newdata = test_df, type = "response") # expected count
# --- Metrics on TEST ----
y <- test_df$travels
# Poisson deviance for predictions μ (handles y=0 safely)
pois_dev <- function(y, mu) {
mu <- pmax(mu, .Machine$[Link])
term <- ifelse(y == 0, -mu, 2 * (y * log(y / mu) - (y - mu)))
# For y>0 formula above gives 2*(y*log(y/mu) - (y-mu));
# For y=0, conventional deviance contribution is 2*mu, but the simplified ifelse yields -mu above.
# Use the standard definition explicitly:
term <- ifelse(y == 0, 2 * mu, 2 * (y * log(y / mu) - (y - mu)))
sum(term)
}
#MSE
sum((y - pred_pois)^2)
sum((y - pred_nb)^2)
sum((y - pred_zip)^2)
#MAE
sum(abs(y - pred_pois))
sum(abs(y - pred_nb))
sum(abs(y - pred_zip))
#Dev_Pois
Dev_Pois = c(pois_dev(y, pred_pois), pois_dev(y, pred_nb), pois_dev(y, pred_zip))
## create the plot for Poisson
ggplot(test_df, aes(x = disc_ons, y = pred_pois, colour = lo_card)) +
geom_point(aes(y = travels), alpha=.5, position=position_jitter(h=.2)) +
geom_line(size = 1) +
labs(x = "Discount During on Season", y = "Expected number of Travels")
## create the plot for NB
ggplot(test_df, aes(x = disc_ons, y = pred_nb, colour = lo_card)) +
geom_point(aes(y = travels), alpha=.5, position=position_jitter(h=.2)) +
geom_line(size = 1) +
labs(x = "Discount During on Season", y = "Expected number of Travels")
## create the plot for ZIP
ggplot(test_df, aes(x = disc_ons, y = pred_zip, colour = lo_card)) +
geom_point(aes(y = travels), alpha=.5, position=position_jitter(h=.2)) +
geom_line(size = 1) +
labs(x = "Discount During on Season", y = "Expected number of Travels")