Assignment 6
#Problem 1
library(readxl)
library(base)
library(stats)
library(dplyr)## Warning: package 'dplyr' was built under R version 4.2.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# Creating a data frame with the provided data
data <- data.frame(
League = c("NFL", "NFL", "NFL", "NFL", "NFL", "AFL", "AFL", "AFL", "AFL", "AFL"),
Distance = c(14.5, 24.5, 34.5, 44.5, 52.0, 14.5, 24.5, 34.5, 44.5, 52.0),
Success = c(68, 74, 61, 38, 2, 62, 49, 43, 25, 7),
Attempts = c(77, 95, 113, 138, 38, 67, 70, 79, 82, 24),
Z = c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
)
# Print the data frame
print(data)## League Distance Success Attempts Z
## 1 NFL 14.5 68 77 0
## 2 NFL 24.5 74 95 0
## 3 NFL 34.5 61 113 0
## 4 NFL 44.5 38 138 0
## 5 NFL 52.0 2 38 0
## 6 AFL 14.5 62 67 1
## 7 AFL 24.5 49 70 1
## 8 AFL 34.5 43 79 1
## 9 AFL 44.5 25 82 1
## 10 AFL 52.0 7 24 1data <- data %>%mutate(Probability=Success/Attempts, Y= ifelse(Probability > 0.5,1,0), Distancesq=Distance^2)
data## League Distance Success Attempts Z Probability Y Distancesq
## 1 NFL 14.5 68 77 0 0.88311688 1 210.25
## 2 NFL 24.5 74 95 0 0.77894737 1 600.25
## 3 NFL 34.5 61 113 0 0.53982301 1 1190.25
## 4 NFL 44.5 38 138 0 0.27536232 0 1980.25
## 5 NFL 52.0 2 38 0 0.05263158 0 2704.00
## 6 AFL 14.5 62 67 1 0.92537313 1 210.25
## 7 AFL 24.5 49 70 1 0.70000000 1 600.25
## 8 AFL 34.5 43 79 1 0.54430380 1 1190.25
## 9 AFL 44.5 25 82 1 0.30487805 0 1980.25
## 10 AFL 52.0 7 24 1 0.29166667 0 2704.00mean(data$Probability)## [1] 0.5296103#(a)
#fitting the logistic model for NFL:
model_NFL <- glm(Y ~ Distance + Distancesq, data = data%>%filter(League== "NFL"),family = binomial)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredmodel_NFL##
## Call: glm(formula = Y ~ Distance + Distancesq, family = binomial, data = data %>%
## filter(League == "NFL"))
##
## Coefficients:
## (Intercept) Distance Distancesq
## -49.0820 7.2407 -0.1494
##
## Degrees of Freedom: 4 Total (i.e. Null); 2 Residual
## Null Deviance: 6.73
## Residual Deviance: 5.416e-10 AIC: 6###ANSWER: The output of the fitted model =>
#log(pi/(1-pi)) = -49.08 + 7.241X - 0.1494X^2
#fitting the logistic model for AFL:
model_AFL = glm(Y~Distance+Distancesq,data=data%>%filter(League == "AFL"), family=binomial)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(model_AFL)##
## Call:
## glm(formula = Y ~ Distance + Distancesq, family = binomial, data = data %>%
## filter(League == "AFL"))
##
## Deviance Residuals:
## 1 2 3 4 5
## 6.788e-06 2.110e-08 1.518e-05 -1.628e-05 -2.110e-08
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.908e+01 3.814e+05 0 1
## Distance 7.241e+00 2.507e+04 0 1
## Distancesq -1.494e-01 3.937e+02 0 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 6.7301e+00 on 4 degrees of freedom
## Residual deviance: 5.4159e-10 on 2 degrees of freedom
## AIC: 6
##
## Number of Fisher Scoring iterations: 23###ANSWER: The output of the fitted model =>
#log(pi/(1-pi)) = -49.08 + 7.241X - 0.1494X^2
#(b)
#fitting the logistic model for AFL:
model = glm(Y~Distance+Distancesq + Z,data=data,family=binomial)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(model)##
## Call:
## glm(formula = Y ~ Distance + Distancesq + Z, family = binomial,
## data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -9.874e-06 -2.110e-08 2.110e-08 4.117e-06 9.208e-06
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.142e+01 4.467e+05 0 1
## Distance 7.567e+00 2.923e+04 0 1
## Distancesq -1.561e-01 4.590e+02 0 1
## Z 1.008e-11 8.594e+04 0 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1.3460e+01 on 9 degrees of freedom
## Residual deviance: 3.9849e-10 on 6 degrees of freedom
## AIC: 8
##
## Number of Fisher Scoring iterations: 24###ANSWER: The output of the fitted model =>
#log(pi/(1-pi)) = -51.42 + 7.567X - 0.156X^2 + 6.63 x 10^-11Z
#(c) Does the quadratic term contribute significantly to the model? The model summary from the above states that p-value is 1 which is greater than 0.05, so quadratic predictor is insignificant
#d) Are the probabilities of scoring field goals from a given distance the same for each league?
#Finding pi equation (as given in the question) at x= 15:
NFL_x = -49.08 + 7.241 *15 - (0.1494* 15^2)
pi_NFL = exp(NFL_x)/(1+exp(NFL_x))
cat("NFL_pi = ", pi_NFL, "\n")## NFL_pi = 1#Finding pi equation (as given in the question) at x = 15:
AFL_x = -49.08 + 7.241 *15 - (0.1494* 15^2)
pi_AFL = exp(AFL_x)/(1+exp(AFL_x))
cat("AFL_pi = ", pi_AFL, "\n")## AFL_pi = 1#Problem 2
library(nnet) # for multinom
library(MASS) # for polr (ordinal regression)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(readxl)
library(dplyr)
path_to_file_12_6 <- "/Users/priyadarshinirajendran/Downloads/table 12.6.xlsx"
path_to_file_12_7 <- "/Users/priyadarshinirajendran/Downloads/tabel 12.7.xlsx"
# Read the Excel files
data_12_6 <- read_excel(path_to_file_12_6)
data_12_7 <- read_excel(path_to_file_12_7)
# Combine the data frames
data_2 <- bind_rows(data_12_6, data_12_7)
print(data_2)## # A tibble: 200 × 5
## Patient RW IR SSPG CC
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0.81 124 55 3
## 2 2 0.95 117 76 3
## 3 3 0.94 143 105 3
## 4 4 1.04 199 108 3
## 5 5 1 240 143 3
## 6 6 0.76 157 165 3
## 7 7 0.91 221 119 3
## 8 8 1.1 186 105 3
## 9 9 0.99 142 98 3
## 10 10 0.78 131 94 3
## # ℹ 190 more rows# Part (a): Multinomial Logistic Regression
# Model without RW
model_multinom_no_RW <- multinom(CC ~ IR + SSPG, data = data_2)## # weights: 12 (6 variable)
## initial value 219.722458
## iter 10 value 104.008931
## iter 20 value 103.538576
## final value 103.538532
## convergedsummary(model_multinom_no_RW)## Call:
## multinom(formula = CC ~ IR + SSPG, data = data_2)
##
## Coefficients:
## (Intercept) IR SSPG
## 2 2.380106 0.01646186 -0.02311931
## 3 6.445230 0.01402701 -0.04358092
##
## Std. Errors:
## (Intercept) IR SSPG
## 2 1.180021 0.003302836 0.005354224
## 3 1.232450 0.003616663 0.006451030
##
## Residual Deviance: 207.0771
## AIC: 219.0771pred_no_RW <- predict(model_multinom_no_RW, data_2)
conf_matrix_no_RW <- table(data_2$CC, pred_no_RW)
classification_rate_no_RW <- sum(diag(conf_matrix_no_RW)) / nrow(data_2)
# Model with RW
model_multinom_RW <- multinom(CC ~ RW + IR + SSPG, data = data_2)## # weights: 15 (8 variable)
## initial value 219.722458
## iter 10 value 100.465877
## iter 20 value 98.755497
## iter 30 value 98.710659
## final value 98.709377
## convergedsummary(model_multinom_RW)## Call:
## multinom(formula = CC ~ RW + IR + SSPG, data = data_2)
##
## Coefficients:
## (Intercept) RW IR SSPG
## 2 -4.922428 8.571885 0.01610192 -0.02946041
## 3 1.021778 6.418457 0.01354776 -0.04804984
##
## Std. Errors:
## (Intercept) RW IR SSPG
## 2 1.266728 1.162236 0.003531896 0.006226657
## 3 1.105563 1.247886 0.003804663 0.007335708
##
## Residual Deviance: 197.4188
## AIC: 213.4188pred_RW <- predict(model_multinom_RW, data_2)
conf_matrix_RW <- table(data_2$CC, pred_RW)
classification_rate_RW <- sum(diag(conf_matrix_RW)) / nrow(data_2)# Part (b): Ordinal Logistic Regression
# Model without RW
library(MASS) # for polr
library(dplyr) # for data manipulation
# Ensure the response variable 'CC' is a factor
data_2$CC <- as.factor(data_2$CC)
# Fit an ordinal logistic regression model without RW
model_ordinal_no_RW <- polr(CC ~ IR + SSPG, data = data_2, method = "logistic")
summary(model_ordinal_no_RW)##
## Re-fitting to get Hessian## Call:
## polr(formula = CC ~ IR + SSPG, data = data_2, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## IR 0.005648 0.001432 3.943
## SSPG -0.028515 0.003033 -9.403
##
## Intercepts:
## Value Std. Error t value
## 1|2 -6.2269 0.6932 -8.9822
## 2|3 -3.4030 0.5278 -6.4480
##
## Residual Deviance: 236.5749
## AIC: 244.5749# Fit an ordinal logistic regression model with RW
model_ordinal_RW <- polr(CC ~ RW + IR + SSPG, data = data_2, method = "logistic")
summary(model_ordinal_RW)##
## Re-fitting to get Hessian## Call:
## polr(formula = CC ~ RW + IR + SSPG, data = data_2, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## RW 2.952429 1.617935 1.825
## IR 0.005207 0.001477 3.526
## SSPG -0.030387 0.003314 -9.170
##
## Intercepts:
## Value Std. Error t value
## 1|2 -3.7376 1.4924 -2.5045
## 2|3 -0.8875 1.4488 -0.6126
##
## Residual Deviance: 233.1457
## AIC: 243.1457model_ordinal_no_RW <- polr(CC ~ IR + SSPG, data = data_2, method = "logistic")
summary(model_ordinal_no_RW)##
## Re-fitting to get Hessian## Call:
## polr(formula = CC ~ IR + SSPG, data = data_2, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## IR 0.005648 0.001432 3.943
## SSPG -0.028515 0.003033 -9.403
##
## Intercepts:
## Value Std. Error t value
## 1|2 -6.2269 0.6932 -8.9822
## 2|3 -3.4030 0.5278 -6.4480
##
## Residual Deviance: 236.5749
## AIC: 244.5749# Model with RW
model_ordinal_RW <- polr(CC ~ RW + IR + SSPG, data = data_2, method = "logistic")
summary(model_ordinal_RW)##
## Re-fitting to get Hessian## Call:
## polr(formula = CC ~ RW + IR + SSPG, data = data_2, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## RW 2.952429 1.617935 1.825
## IR 0.005207 0.001477 3.526
## SSPG -0.030387 0.003314 -9.170
##
## Intercepts:
## Value Std. Error t value
## 1|2 -3.7376 1.4924 -2.5045
## 2|3 -0.8875 1.4488 -0.6126
##
## Residual Deviance: 233.1457
## AIC: 243.1457#Problem 3
library(ggplot2)
library(MASS)
library(dplyr)# Create the data frame
data <- data.frame(Y = c(11, 7, 7, 19, 9, 4, 3, 1, 3),
N = c(0.0950, 0.1920, 0.0750, 0.2078, 0.1382, 0.0540, 0.1292, 0.0503, 0.0629))
# Fit the models
model_ls <- lm(Y ~ N, data = data)
model_tls <- lm(log(Y + 1) ~ N, data = data)
model_poisson <- glm(Y ~ N, family = poisson, data = data)
# Calculate residuals for each model
data <- data %>%
mutate(resid_ls = residuals(model_ls),
resid_tls = residuals(model_tls),
resid_poisson = residuals(model_poisson, type = "response"))
# Plotting residuals
ggplot(data, aes(x = N)) +
geom_point(aes(y = resid_ls), color = "blue") +
geom_point(aes(y = resid_tls), color = "green") +
geom_point(aes(y = resid_poisson), color = "red") +
labs(title = "Residuals Comparison", x = "N", y = "Residuals") +
theme_minimal() +
geom_hline(yintercept = 0, linetype = "dashed")
###Based on the above image Transformed Least Square gives the best
description of the given data
#Problem 4
library(stats)
library(MASS)
# Data
magazines <- c('Cosmopolitan', 'Redbook', 'Glamour', 'Southern Living', 'Vogue',
'Sunset', 'House and Garden', 'New York Magazine', 'House Beautiful',
'Mademoiselle', 'Psychology Today', 'Life Magazine', 'Smithsonian',
'Rolling Stone', 'Modern Bride', 'Parents', 'Architectural Digest',
'Harper Bazaar', 'Apartment Life', 'Bon Appetit', 'Gourmet',
'Town and Country', 'True Story', 'Brides', 'Book Digest Magazine',
'W', 'Yankee', 'Playgirl', 'Saturday Review', 'New Woman', 'Ms.',
'Cuisine', 'Mother Earth News', '1001 Decorating Ideas', 'Self',
'Decorating Craft Ideas', 'Saturday Evening Post', 'McCall Needlework and Craft',
'Weight Watchers', 'High Times', 'Soap Opera Digest')
P <- c(25, 15, 20, 17, 23, 17, 14, 22, 12, 15, 8, 7, 9, 12, 1, 6, 12, 9, 7, 9, 7, 1, 77, 13, 5, 7, 13, 4, 6, 3, 6, 4, 3, 3, 5, 4, 4, 3, 3, 4, 2)
R <- c(50, 49.7, 34, 30.7, 27, 26.3, 24.6, 16.9, 16.7, 14.6, 13.8, 13.2, 13.1, 10.6, 8.8, 8.7, 8.5, 8.3, 8.2, 8.2, 7.3, 7, 6.6, 6.2, 5.8, 5.1, 4.1, 3.9, 3.9, 3.5, 3.3, 3, 2.5, 2.3, 2.3, 1.8, 1.5, 1.3, 1.3, 1, 0.3)
# Combine into data frame
data <- data.frame(Magazine = magazines, P = P, R = R)
# Log transform
data$log_P <- log(data$P)
data$log_R <- log(data$R)
# Fit least squares regression
ls_fit <- lm(log_R ~ log_P, data = data)
# Identify special features of observations 15, 22, 23, and 41
problematic_indices <- c(15, 22, 23, 41)
problematic_points <- data[problematic_indices,]
# Fit robust regression
robust_fit <- rlm(log_R ~ log_P, data = data)
# Compare robust fit to least squares fit after removing problematic observations
data_reduced <- data[-problematic_indices,]
ls_fit_reduced <- lm(log_R ~ log_P, data = data_reduced)
# Output the summaries for comparison
summary(ls_fit)##
## Call:
## lm(formula = log_R ~ log_P, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.01534 -0.53524 0.04836 0.50718 1.94245
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2323 0.3398 0.684 0.498
## log_P 0.8354 0.1571 5.318 4.57e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8719 on 39 degrees of freedom
## Multiple R-squared: 0.4203, Adjusted R-squared: 0.4055
## F-statistic: 28.28 on 1 and 39 DF, p-value: 4.575e-06summary(robust_fit)##
## Call: rlm(formula = log_R ~ log_P, data = data)
## Residuals:
## Min 1Q Median 3Q Max
## -2.612335 -0.475388 -0.002696 0.453661 2.492415
##
## Coefficients:
## Value Std. Error t value
## (Intercept) -0.3177 0.2755 -1.1532
## log_P 1.1090 0.1273 8.7083
##
## Residual standard error: 0.7048 on 39 degrees of freedomsummary(ls_fit_reduced)##
## Call:
## lm(formula = log_R ~ log_P, data = data_reduced)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2759 -0.3202 0.1032 0.3523 1.0137
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9941 0.2845 -3.494 0.00131 **
## log_P 1.4351 0.1318 10.891 8.67e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5185 on 35 degrees of freedom
## Multiple R-squared: 0.7722, Adjusted R-squared: 0.7657
## F-statistic: 118.6 on 1 and 35 DF, p-value: 8.667e-13