#Coded by Tawkir Ahmed
library(ggplot2) # USed for plotting data
library(dplyr) # Used to extract columns in the data
library(rms) # Used to extract p-value from logistic model
library(aod)
theme_set(theme_gray() ) # the default
# logistic model
covid <- read.csv(choose.files())
labs <- attributes(covid)$labels
summary(covid)
# collapse all missing values to NA
covid$x <- factor(covid$x, levels=c("1", "2", "3", "4", "5", "6", "7"))
# run our regression model
before <- glm(x~ x1+ x2+ x3+ x4+
x5+ x6,
data=covid, family="binomial")
summary(before)
# check that model is good fit or not
with(before, cbind(res.deviance = deviance, df = df.residual,
p = pchisq(deviance, df.residual, lower.tail=FALSE)))
## odds ratios and 95% CI ***********
exp(cbind(OR = coef(before), confint(before)))
# MORE SUMMARIES ###########################################
anova(before) # Coefficients w/inferential tests
coef(before) # Coefficients
hist(coef(before))
confint(before) # CI for coefficients
hist(confint(before))
resid(before) # Residuals case-by-case
hist(residuals(before),main="Before COVID 19" ) # Histogram of residuals
plot(residuals(before), main="Before COVID 19" )
logLik(before)
BIC(before)
PseudoR2(before)
predict(before)
hist(predict(before))
#peseudo r square
before$null.deviance
before$deviance
modelChi <- before$null.deviance - before$deviance
pseudo <- modelChi / before$null.deviance
pseudo
# Compute the pseudo p-value
Chidf <- before$df.null - before$df.residual
modelChi <- before$null.deviance - before$deviance
1 - pchisq(modelChi, Chidf)
#RSS(residual sum of square)
RSS <- c(crossprod(before$residuals))
RSS
#Mean square error
MSE <- RSS / length(before$residuals)
MSE
#Root MSE
RMSE <- sqrt(MSE)
RMSE
#Pearson estimated residual variance
sig2 <- RSS / before$df.residual
sig2
# CLEAN UP #################################################
plot(before) #plot 4 times to get four graph
library(effects)
x1 <- predictorEffect("x6", before)
plot(x1)
x2 <- predictorEffect("x2", before)
plot(x2, main="Before COVID 19",xlab="X2", ylab="x") #single plot
plot(allEffects(before)) #all plot
0 comments:
Post a Comment