Exam1.1 {StroupGLMM}

Example1.1 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup(p-5)

Description

Exam1.1 is used for inspecting probability distribution and to define a plausible process through linear models and generalized linear models.

Author(s)

1. Muhammad Yaseen ([e-mail])

2. Adeela Munawar ([e-mail])

References

1. Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.

See Also

Table1.1

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
    # , subset
    # , weights
    # , na.action
    , method      =  "qr"
    , model       =  TRUE
    , x           =  FALSE
    , y           =  FALSE
    , qr          =  TRUE
    , singular.ok =  TRUE
    , contrasts   =  NULL
    # , offset
    # , ...
  )
summary(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
    , weights     =  NULL
    # , subset
    # , na.action
    , start       =  NULL
    # , etastart
    # , mustart
    # , offset
    # , control     =  list(...)
    # , model       =  TRUE
    , method      =  "glm.fit"
    , x           =  FALSE
    , y           =  TRUE
    , contrasts   =  NULL
    # , ...
  )
## this glm() function gives warning message of non-integer success
summary(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
    , weights   =  NULL
    # , subset
    # , na.action
    , start     =  NULL
    # , etastart
    # , mustart
    # , offset
    # , control   =  list(...)
    # , model     =  TRUE
    , method    =  "glm.fit"
    , x         =  FALSE
    , y         =  TRUE
    , contrasts =  NULL
    # , ...
  )
## problem of "warning message of non-integer success" is overome by using quasibinomial family
summary(Exam1.1.glm2)

#-------------------------------------------------------------
## GLM fitting with survey package(produces same result as using quasi binomial family in glm)
#-------------------------------------------------------------
library(survey)
design   <-
  svydesign(
      ids          =  ~1
    , probs        =  NULL
    , strata       =  NULL
    , variables    =  NULL
    , fpc          =  NULL
    , data         =  Table1.1
    # , nest         =  FALSE
    # , check.strata =  !nest
    , weights      =  NULL
    , pps          =  FALSE
    # , ...
  )

Exam1.1.svyglm  <-
  svyglm(
      formula  =  y/Nx~x
    , design   =  design
    # , ...
    , family   =  quasibinomial(link="logit")
  )
# summary(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)
)