R/Exam1.1.R

In myaseen208/StroupGLMM: R Codes and Datasets for Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup

#' @title Example1.1 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup(p-5)
#' @name   Exam1.1
#' @description Exam1.1 is used for inspecting probability distribution and to define a plausible process through
#' linear models and generalized linear models.
#'
#' @author
#'  \enumerate{
#'          \item  Muhammad Yaseen (\email{myaseen208@gmail.com})
#'          \item Adeela Munawar (\email{adeela.uaf@gmail.com})
#'          }
#'
#' @references
#' \enumerate{
#'    \item Stroup, W. W. (2012).
#'      \emph{Generalized Linear Mixed Models: Modern Concepts, Methods and Applications}.
#'        CRC Press.
#'          }
#'
#' @seealso
#'    \link{Table1.1}
#'
#' @import parameters
#' @importFrom ggplot2 ggplot
#' @importFrom stats lm summary.lm glm summary.glm cor
#' @importFrom survey svydesign svyglm
#'
#' @examples
#' #-------------------------------------------------------------
#' ## Linear Model and results discussed in Article 1.2.1 after Table1.1
#' #-------------------------------------------------------------
#' data(Table1.1)
#' Exam1.1.lm1 <- lm(formula =  y/Nx ~ x, data = Table1.1)
#' summary(Exam1.1.lm1 )
#' library(parameters)
#' model_parameters(Exam1.1.lm1)
#'
#' #-------------------------------------------------------------
#' ## GLM fitting with logit link (family=binomial)
#' #-------------------------------------------------------------
#' Exam1.1.glm1 <-
#'         glm(
#'               formula = y/Nx ~ x
#'             , family  =  binomial(link = "logit")
#'             , data    = Table1.1
#'             )
#' ## this glm() function gives warning message of non-integer success
#' summary(Exam1.1.glm1)
#' model_parameters(Exam1.1.glm1)
#'
#' #-------------------------------------------------------------
#' ## GLM fitting with logit link (family = Quasibinomial)
#' #-------------------------------------------------------------
#' Exam1.1.glm2 <-
#'       glm(
#'            formula = y/Nx~x
#'          , family =  quasibinomial(link = "logit")
#'          , data   =  Table1.1
#'          )
#' ## problem of "warning message of non-integer success" is overome by using quasibinomial family
#' summary(Exam1.1.glm2)
#' model_parameters(Exam1.1.glm2)
#'
#' #-------------------------------------------------------------
#' ## GLM fitting with survey package(produces same result as using quasi binomial family in glm)
#' #-------------------------------------------------------------
#' library(survey)
#' design   <- svydesign(ids =  ~1, data =  Table1.1)
#'
#' Exam1.1.svyglm  <-
#'         svyglm(
#'                  formula  =  y/Nx~x
#'                , design   =  design
#'                , family   =  quasibinomial(link = "logit")
#'                )
#' summary(Exam1.1.svyglm)
#' model_parameters(Exam1.1.svyglm)
#'
#' #-------------------------------------------------------------
#' ## Figure 1.1
#' #-------------------------------------------------------------
#' Newdata     <-
#'   data.frame(
#'     Table1.1
#'     , LM       =  Exam1.1.lm1$fitted.values
#'     , GLM      =  Exam1.1.glm1$fitted.values
#'     , QB       =  Exam1.1.glm2$fitted.values
#'     , SM       =  Exam1.1.svyglm$fitted.values
#'   )
#' #-------------------------------------------------------------
#' ## One Method to plot  Figure1.1
#' #-------------------------------------------------------------
#' library(ggplot2)
#'
#' Figure1.1   <-
#'   ggplot(
#'       data     = Newdata
#'     , mapping  = aes(x = x, y = y/Nx)
#'   )     +
#'   geom_point (
#'     mapping  = aes(colour = "black")
#'   )  +
#'   geom_point (
#'     data     = Newdata
#'     , mapping  = aes(x = x, y = LM, colour = "blue"), shape = 2
#'   )  +
#'   geom_line(
#'     data     = Newdata
#'     , mapping  = aes(x = x, y = LM, colour = "blue")
#'   )   +
#'   geom_point (
#'     data     = Newdata
#'     , mapping  = aes(x = x, y = GLM, colour ="red"), shape = 3
#'   ) +
#'   geom_smooth (
#'     data     = Newdata
#'     , mapping  = aes(x = x, y = GLM, colour = "red")
#'     , stat     = "smooth"
#'   ) +
#'   theme_bw()    +
#'   scale_colour_manual (
#'     values = c("black", "blue", "red"),
#'     labels = c("observed", "LM", "GLM")
#'   )  +
#'   guides (
#'     colour   = guide_legend(title = "Plot")
#'   ) +
#'   labs (
#'     title     = "Linear Model vs Logistic Model"
#'   ) +
#'   labs (
#'     y         = "p"
#'   )
#' print(Figure1.1)
#'
#' #-------------------------------------------------------------
#' ## Another way to plot Figure 1.1
#' #-------------------------------------------------------------
#' newdata   <-
#'   data.frame(
#'     P     =  c(
#'                 Table1.1$y/Table1.1$Nx
#'               , Exam1.1.lm1$fitted.values
#'               , Exam1.1.glm1$fitted.values
#'                )
#'     , X     =  rep(Table1.1$x, 3)
#'     , group =  rep(c('Obs','LM','GLM'), each = length(Table1.1$x))
#'   )
#'
#' Figure1.1      <-
#'   ggplot(
#'       data    = newdata
#'     , mapping = aes(x = X , y = P)
#'   )    +
#'   geom_point(
#'     mapping = aes(x = X , y = P, colour = group , shape=group)
#'   ) +
#'   geom_smooth(
#'     data    = subset(x = newdata, group == "LM")
#'     , mapping = aes(x=X,y=P)
#'     , col     = "green"
#'   ) +
#'   geom_smooth(
#'     data    = subset(x = newdata, group=="GLM")
#'     , mapping = aes(x = X , y = P)
#'     , col     = "red"
#'   ) +
#'   theme_bw() +
#'   labs(
#'     title   = "Linear Model vs Logistic Model"
#'   )
#' print(Figure1.1)
#'
#' #-------------------------------------------------------------
#' ## Correlation among p and fitted values using Gaussian link
#' #-------------------------------------------------------------
#' (lmCor <- cor(Table1.1$y/Table1.1$Nx, Exam1.1.lm1$fitted.values))
#'
#' #-------------------------------------------------------------
#' ## Correlation among p and fitted values using quasi binomial link
#' #-------------------------------------------------------------
#' (glmCor <- cor(Table1.1$y/Table1.1$Nx, Exam1.1.glm1$fitted.values))
NULL