Friday, November 13, 2020

Logistic Model

 #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