Problem 1
# manually input the data from table 13.2
data <- data.frame(
Household = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20),
Income = c(38000, 51200, 39600, 43400, 47700, 53000, 41500, 40800, 45400, 52400, 38700, 40100, 49500, 38000, 42000, 54000, 51700, 39400, 40900, 52800),
Ownership = c(0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1))
# set seed
set.seed(123)
# randomly sample 15 rows and all columns
sample_data <- data[sample(1:nrow(data), 15), ]
# view sample
print(sample_data)
## Household Income Ownership
## 15 15 42000 1
## 19 19 40900 0
## 14 14 38000 0
## 3 3 39600 0
## 10 10 52400 1
## 2 2 51200 1
## 6 6 53000 0
## 11 11 38700 1
## 5 5 47700 0
## 4 4 43400 1
## 18 18 39400 0
## 9 9 45400 1
## 20 20 52800 1
## 17 17 51700 1
## 12 12 40100 0
A. Logistic Regression
# fit logistic regression model
model1 <- glm(Ownership ~ Income, data = sample_data, family = binomial)
summary(model1)
##
## Call:
## glm(formula = Ownership ~ Income, family = binomial, data = sample_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.9641123 4.7647986 -1.462 0.144
## Income 0.0001583 0.0001064 1.488 0.137
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 20.728 on 14 degrees of freedom
## Residual deviance: 18.176 on 13 degrees of freedom
## AIC: 22.176
##
## Number of Fisher Scoring iterations: 4
# report estimated standard errors
se_model1 <- summary(model1)$coefficients[, "Std. Error"]
cat("Standard Errors - Linear Model:\n")
## Standard Errors - Linear Model:
print(se_model1)
## (Intercept) Income
## 4.7647986325 0.0001063738
B. Assess Model Deviance
deviance_model1 <- deviance(model1)
df_resid_model1 <- df.residual(model1)
p_value <- 1 - pchisq(deviance_model1, df_resid_model1)
cat("Model deviance:", deviance_model1, "\n")
## Model deviance: 18.17622
cat("Degrees of freedom:", df_resid_model1, "\n")
## Degrees of freedom: 13
cat("p-value:", p_value, "\n")
## p-value: 0.1509425
C.Beta Interpretation
beta_1 <- coef(model1)[2]
se_beta_1 <- summary(model1)$coefficients[2,2]
odds_ratio <- exp(beta_1)
cat("β1 =", beta_1, "with SE =", se_beta_1, "\n")
## β1 = 0.0001583349 with SE = 0.0001063738
cat("Odds ratio for a $1000 increase in income =", exp(1000 * beta_1), "\n")
## Odds ratio for a $1000 increase in income = 1.171559
D. Quadratic Term
# add quadratic term and test for significance
sample_data$Income_sq <- sample_data$Income^2
model2 <- glm(Ownership ~ Income + Income_sq, data = sample_data, family = binomial)
summary(model2)
##
## Call:
## glm(formula = Ownership ~ Income + Income_sq, family = binomial,
## data = sample_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.304e+01 6.906e+01 -0.913 0.361
## Income 2.641e-03 3.043e-03 0.868 0.386
## Income_sq -2.707e-08 3.307e-08 -0.819 0.413
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 20.728 on 14 degrees of freedom
## Residual deviance: 17.454 on 12 degrees of freedom
## AIC: 23.454
##
## Number of Fisher Scoring iterations: 4
# report estimated standard errors
se_model2 <- summary(model2)$coefficients[, "Std. Error"]
cat("\nStandard Errors - Quadratic Model:\n")
##
## Standard Errors - Quadratic Model:
print(se_model2)
## (Intercept) Income Income_sq
## 6.906031e+01 3.042840e-03 3.306579e-08
# likelihood ratio test to determine if quadratic term improves fit
anova_result <- anova(model1, model2, test = "Chisq")
cat("\nLikelihood Ratio Test between Linear and Quadratic Model:\n")
##
## Likelihood Ratio Test between Linear and Quadratic Model:
print(anova_result)
## Analysis of Deviance Table
##
## Model 1: Ownership ~ Income
## Model 2: Ownership ~ Income + Income_sq
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 13 18.176
## 2 12 17.454 1 0.72221 0.3954
Problem 2
# manually input the data from table 13.4
data2 <- data.frame(
Discount = c(5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25),
SampleSize = c(500, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500),
Redeemed = c(100, 122, 147, 176, 211, 244, 277, 310, 343, 372, 391))
# set seed
set.seed(123)
# randomly sample 15 rows and all columns
sample_data2 <- data2[sample(1:nrow(data2), 8), ]
# view sample
print(sample_data2)
## Discount SampleSize Redeemed
## 3 9 500 147
## 11 25 500 391
## 2 7 500 122
## 6 15 500 244
## 10 23 500 372
## 5 13 500 211
## 4 11 500 176
## 9 21 500 343
A. Logistic Regression
# compute proportion redeemed
sample_data2$Prop <- sample_data2$Redeemed / sample_data2$SampleSize
# Fit logistic regression model (logit link, linear predictor)
model3 <- glm(cbind(Redeemed, SampleSize - Redeemed) ~ Discount,
family = binomial, data = sample_data2)
summary(model1)
##
## Call:
## glm(formula = Ownership ~ Income, family = binomial, data = sample_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.9641123 4.7647986 -1.462 0.144
## Income 0.0001583 0.0001064 1.488 0.137
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 20.728 on 14 degrees of freedom
## Residual deviance: 18.176 on 13 degrees of freedom
## AIC: 22.176
##
## Number of Fisher Scoring iterations: 4
# report estimated standard errors
se_model3 <- summary(model3)$coefficients[, "Std. Error"]
cat("Standard Errors - Linear Model:\n")
## Standard Errors - Linear Model:
print(se_model3)
## (Intercept) Discount
## 0.09399535 0.00574384
B. Assess Model Deviance & Adequacy
deviance_model3 <- deviance(model3)
df_resid_model3 <- df.residual(model3)
p_value_model3 <- 1 - pchisq(deviance_model3, df_resid_model3)
cat("Model Deviance:", deviance_model3, "\n")
## Model Deviance: 0.2439253
cat("Residual df:", df_resid_model3, "\n")
## Residual df: 6
cat("p-value:", p_value_model3, "\n")
## p-value: 0.999724
C.Graph Data and Model
ggplot(sample_data2, aes(x = Discount, y = Prop)) +
geom_point(size = 3) +
stat_smooth(method = "glm",
method.args = list(family = "binomial"),
se = FALSE, color = "blue", formula = y ~ x) +
ggtitle("Logistic Regression Fit (Linear Predictor)") +
ylab("Proportion Redeemed") +
theme_minimal()
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
D. Quadratic Term
# add quadratic term and test for significance
sample_data2$Discount_sq <- sample_data2$Discount^2
model4 <- glm(cbind(Redeemed, SampleSize - Redeemed) ~ Discount + Discount_sq,
family = binomial, data = sample_data2)
summary(model4)
##
## Call:
## glm(formula = cbind(Redeemed, SampleSize - Redeemed) ~ Discount +
## Discount_sq, family = binomial, data = sample_data2)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1318923 0.2939898 -7.252 4.12e-13 ***
## Discount 0.1416884 0.0404542 3.502 0.000461 ***
## Discount_sq -0.0001645 0.0012372 -0.133 0.894254
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 645.87550 on 7 degrees of freedom
## Residual deviance: 0.22625 on 5 degrees of freedom
## AIC: 58.155
##
## Number of Fisher Scoring iterations: 3
# report estimated standard errors
se_model4 <- summary(model4)$coefficients[, "Std. Error"]
cat("\nStandard Errors - Quadratic Model:\n")
##
## Standard Errors - Quadratic Model:
print(se_model4)
## (Intercept) Discount Discount_sq
## 0.293989841 0.040454177 0.001237214
# likelihood ratio test to determine if quadratic term improves fit
anova_result2 <- anova(model3, model4, test = "Chisq")
cat("\nLikelihood Ratio Test between Linear and Quadratic Model:\n")
##
## Likelihood Ratio Test between Linear and Quadratic Model:
print(anova_result2)
## Analysis of Deviance Table
##
## Model 1: cbind(Redeemed, SampleSize - Redeemed) ~ Discount
## Model 2: cbind(Redeemed, SampleSize - Redeemed) ~ Discount + Discount_sq
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 6 0.24393
## 2 5 0.22625 1 0.017671 0.8942
E. Plot Both Models
# create prediction grid
discount_seq <- seq(min(sample_data2$Discount), max(sample_data2$Discount), length.out = 100)
# for linear model (only needs Discount)
newdata_linear <- data.frame(Discount = discount_seq)
newdata_linear$LinearFit <- predict(model3, newdata = newdata_linear, type = "response")
# for quadratic model (needs both Discount and Discount_sq)
newdata_quad <- data.frame(
Discount = discount_seq,
Discount_sq = discount_seq^2
)
newdata_quad$QuadraticFit <- predict(model4, newdata = newdata_quad, type = "response")
# combine into one plot
ggplot(sample_data2, aes(x = Discount, y = Prop)) +
geom_point(size = 3, color = "black") +
geom_line(data = newdata_linear, aes(x = Discount, y = LinearFit),
color = "blue", linetype = "dashed", size = 1) +
geom_line(data = newdata_quad, aes(x = Discount, y = QuadraticFit),
color = "red", size = 1) +
ggtitle("Linear vs Quadratic Logistic Regression Fit") +
ylab("Proportion Redeemed") +
theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
F. Wald Statistic
summary_model4 <- summary(model4)
wald_stats <- summary_model4$coefficients[, "Estimate"] / summary_model4$coefficients[, "Std. Error"]
cat("Wald Statistics (Quadratic Model):\n")
## Wald Statistics (Quadratic Model):
print(wald_stats)
## (Intercept) Discount Discount_sq
## -7.2515849 3.5024418 -0.1329235
G. 95% CI
confint_default <- confint.default(model4) # uses normal approx
cat("\n95% CI for Model Parameters (Quadratic Model):\n")
##
## 95% CI for Model Parameters (Quadratic Model):
print(confint_default)
## 2.5 % 97.5 %
## (Intercept) -2.708101801 -1.55568280
## Discount 0.062399669 0.22097713
## Discount_sq -0.002589349 0.00226044
# make sure the quadratic term exists
sample_data2$Discount_sq <- sample_data2$Discount^2
# function to fit the quadratic model on a bootstrap resample and return coefs
boot_fn <- function(d, i) {
dd <- d[i, ] # resample rows with replacement
dd$Discount_sq <- dd$Discount^2 # ensure term present in resample
fit <- glm(cbind(Redeemed, SampleSize - Redeemed) ~ Discount + Discount_sq,
family = binomial, data = dd)
coef(fit)
}
set.seed(123)
b <- boot(data = sample_data2, statistic = boot_fn, R = 2000)
# bootstrap percentile and BCa 95% CIs for each parameter
ci_intercept_perc <- boot.ci(b, type = "perc", index = 1)
ci_discount_perc <- boot.ci(b, type = "perc", index = 2)
ci_quad_perc <- boot.ci(b, type = "perc", index = 3)
ci_intercept_bca <- boot.ci(b, type = "bca", index = 1)
ci_discount_bca <- boot.ci(b, type = "bca", index = 2)
ci_quad_bca <- boot.ci(b, type = "bca", index = 3)
# compare to Wald (normal-approx) CIs from model4
wald_ci <- confint.default(model4)
cat("\nWald 95% CIs (normal approx):\n")
##
## Wald 95% CIs (normal approx):
print(wald_ci)
## 2.5 % 97.5 %
## (Intercept) -2.708101801 -1.55568280
## Discount 0.062399669 0.22097713
## Discount_sq -0.002589349 0.00226044
cat("\nBootstrap 95% CIs (Percentile):\n")
##
## Bootstrap 95% CIs (Percentile):
print(list(
Intercept = ci_intercept_perc$percent[4:5],
Discount = ci_discount_perc$percent[4:5],
Discount2 = ci_quad_perc$percent[4:5]
))
## $Intercept
## [1] -2.395187 -1.996748
##
## $Discount
## [1] 0.1192051 0.1737433
##
## $Discount2
## [1] -0.0010739434 0.0006001724
cat("\nBootstrap 95% CIs (BCa):\n")
##
## Bootstrap 95% CIs (BCa):
print(list(
Intercept = ci_intercept_bca$bca[4:5],
Discount = ci_discount_bca$bca[4:5],
Discount2 = ci_quad_bca$bca[4:5]
))
## $Intercept
## [1] -2.426496 -2.007142
##
## $Discount
## [1] 0.1207236 0.1787630
##
## $Discount2
## [1] -0.0011194439 0.0005952573
Problem 3
# manually input the data from table 13.4
data3 <- data.frame(
Valve = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15),
Failures = c(5, 3, 0, 1, 4, 0, 0, 1, 0, 0, 0, 1, 0, 7, 0),
Months = c(18, 15, 11, 14, 23, 10, 5, 8, 7, 12, 3, 7, 2, 30, 9))
# set seed
set.seed(123)
# randomly sample 15 rows and all columns
sample_data3 <- data3[sample(1:nrow(data3), 6), ]
# view sample
print(sample_data3)
## Valve Failures Months
## 15 15 0 9
## 3 3 0 11
## 14 14 7 30
## 10 10 0 12
## 2 2 3 15
## 6 6 0 10
A. Poisson Model
# fit Poisson regression model
poisson_model <- glm(Failures ~ Months, data = sample_data3, family = poisson)
summary(poisson_model)
##
## Call:
## glm(formula = Failures ~ Months, family = poisson, data = sample_data3)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.42400 1.07492 -2.255 0.024130 *
## Months 0.14791 0.04029 3.671 0.000241 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 23.6179 on 5 degrees of freedom
## Residual deviance: 6.8794 on 4 degrees of freedom
## AIC: 17.679
##
## Number of Fisher Scoring iterations: 6
# extract standard errors
se_poisson <- summary(poisson_model)$coefficients[, "Std. Error"]
cat("Estimated Standard Errors of Coefficients:\n")
## Estimated Standard Errors of Coefficients:
print(se_poisson)
## (Intercept) Months
## 1.07491872 0.04028828
# define bootstrap function for Poisson regression
boot_poisson <- function(data, indices) {
d <- data[indices, ] # resample rows
model <- glm(Failures ~ Months, family = poisson, data = d)
return(coef(model)) # return coefficients: intercept and slope
}
# set seed for reproducibility
set.seed(123)
# perform bootstrap with 1000 replications
boot_results <- boot(data = sample_data3, statistic = boot_poisson, R = 1000)
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
# View basic results
print(boot_results)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = sample_data3, statistic = boot_poisson, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -2.4239972 -38.32675 43.557124
## t2* 0.1479117 2.12340 3.070857
# Compute 95% percentile-based CIs
boot.ci(boot_results, type = "perc", index = 1) # Intercept
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_results, type = "perc", index = 1)
##
## Intervals :
## Level Percentile
## 95% (-121.397, -0.758 )
## Calculations and Intervals on Original Scale
boot.ci(boot_results, type = "perc", index = 2) # Months
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_results, type = "perc", index = 2)
##
## Intervals :
## Level Percentile
## 95% ( 0.0000, 8.1664 )
## Calculations and Intervals on Original Scale
B. Assess Model Deviance & Adequacy
model_dev <- deviance(poisson_model)
df_resid <- df.residual(poisson_model)
p_val <- 1 - pchisq(model_dev, df_resid)
cat("Model Deviance:", model_dev, "\n")
## Model Deviance: 6.879417
cat("Residual Degrees of Freedom:", df_resid, "\n")
## Residual Degrees of Freedom: 4
cat("p-value:", p_val, "\n")
## p-value: 0.1423994
C.Graph Fitted vs Observed Data
# add predicted values to the dataset
sample_data3$Fitted <- predict(poisson_model, type = "response")
# plot observed vs fitted values
ggplot(sample_data3, aes(x = Months)) +
geom_point(aes(y = Failures), size = 3, color = "black") +
geom_line(aes(y = Fitted), color = "blue", size = 1) +
ggtitle("Observed vs Fitted Poisson Regression") +
ylab("Number of Failures") +
xlab("Months in Service") +
theme_minimal()
