Multivariate Covariance Generalized Linear Models

---

To install the stable version of mcglm, use devtools::install_git(). For more information, visit mcglm/README.

library(devtools)
install_git("http://git.leg.ufpr.br/wbonat/mcglm.git")

library(mcglm)
packageVersion("mcglm")

## [1] '0.0.2'

---

The Australian Health Survey

##----------------------------------------------------------------------
## Loadin the Australian Health Survey data.

data(ahs, package="mcglm")

## Object structure.
str(ahs)

## 'data.frame':    5190 obs. of  15 variables:
##  $ sex     : Factor w/ 2 levels "male","female": 1 1 1 1 1 1 1 1 1 1 ...
##  $ age     : num  0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 ...
##  $ income  : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ levyplus: int  0 0 0 0 0 0 0 0 0 1 ...
##  $ freepoor: int  0 0 0 0 0 0 1 1 1 0 ...
##  $ freerepa: int  0 0 0 0 0 1 0 0 0 0 ...
##  $ illness : int  0 1 1 2 3 3 0 1 1 0 ...
##  $ actdays : int  0 4 0 0 0 0 0 0 1 0 ...
##  $ hscore  : int  0 0 1 3 0 3 0 0 3 1 ...
##  $ chcond  : Factor w/ 3 levels "otherwise","not limited",..: 1 1 1 1 1 2 1 2 1 1 ...
##  $ Ndoc    : int  0 3 0 0 0 0 0 1 0 0 ...
##  $ Nndoc   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Nadm    : int  0 0 0 0 0 0 0 1 0 0 ...
##  $ Nhosp   : int  0 0 0 0 0 0 0 5 0 0 ...
##  $ Nmed    : int  1 1 1 1 2 0 0 1 0 1 ...

## Descriptive measures.
summary(ahs)

##      sex            age             income          levyplus     
##  male  :2488   Min.   :0.1900   Min.   :0.0000   Min.   :0.0000  
##  female:2702   1st Qu.:0.2200   1st Qu.:0.2500   1st Qu.:0.0000  
##                Median :0.3200   Median :0.5500   Median :0.0000  
##                Mean   :0.4064   Mean   :0.5832   Mean   :0.4428  
##                3rd Qu.:0.6200   3rd Qu.:0.9000   3rd Qu.:1.0000  
##                Max.   :0.7200   Max.   :1.5000   Max.   :1.0000  
##     freepoor          freerepa         illness     
##  Min.   :0.00000   Min.   :0.0000   Min.   :0.000  
##  1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.000  
##  Median :0.00000   Median :0.0000   Median :1.000  
##  Mean   :0.04277   Mean   :0.2102   Mean   :1.432  
##  3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:2.000  
##  Max.   :1.00000   Max.   :1.0000   Max.   :5.000  
##     actdays            hscore               chcond    
##  Min.   : 0.0000   Min.   : 0.000   otherwise  :2493  
##  1st Qu.: 0.0000   1st Qu.: 0.000   not limited:2092  
##  Median : 0.0000   Median : 0.000   limited    : 605  
##  Mean   : 0.8619   Mean   : 1.218                     
##  3rd Qu.: 0.0000   3rd Qu.: 2.000                     
##  Max.   :14.0000   Max.   :12.000                     
##       Ndoc            Nndoc              Nadm       
##  Min.   :0.0000   Min.   : 0.0000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.: 0.0000   1st Qu.:0.0000  
##  Median :0.0000   Median : 0.0000   Median :0.0000  
##  Mean   :0.3017   Mean   : 0.2146   Mean   :0.1736  
##  3rd Qu.:0.0000   3rd Qu.: 0.0000   3rd Qu.:0.0000  
##  Max.   :9.0000   Max.   :11.0000   Max.   :5.0000  
##      Nhosp             Nmed      
##  Min.   : 0.000   Min.   :0.000  
##  1st Qu.: 0.000   1st Qu.:0.000  
##  Median : 0.000   Median :1.000  
##  Mean   : 1.334   Mean   :1.218  
##  3rd Qu.: 0.000   3rd Qu.:2.000  
##  Max.   :80.000   Max.   :8.000

##----------------------------------------------------------------------
## Frequency tables.

names(ahs)[c(1, 4:7, 10)]

## [1] "sex"      "levyplus" "freepoor" "freerepa" "illness" 
## [6] "chcond"

par(mfrow=c(2,3))
## sapply(ahs[, c(1, 4:7, 10)],
##        FUN=function(x){
##            ## pie(table(x))
##            barplot(prop.table(table(x)))
##        })

barplot(prop.table(xtabs(~sex, data=ahs)),
        ylab="Sample proportion",
        xlab="Sex")

barplot(prop.table(xtabs(~levyplus, data=ahs)),
        ylab="Sample proportion",
        xlab="levyplus")

barplot(prop.table(xtabs(~freepoor, data=ahs)),
        ylab="Sample proportion",
        xlab="freepoor")

barplot(prop.table(xtabs(~freerepa, data=ahs)),
        ylab="Sample proportion",
        xlab="freerepa")

barplot(prop.table(xtabs(~illness, data=ahs)),
        ylab="Sample proportion",
        xlab="illness")

barplot(prop.table(xtabs(~chcond, data=ahs)),
        ylab="Sample proportion",
        xlab="chcond")

layout(1)

xt <- xtabs(~age+sex, data=ahs)
mosaicplot(xt)

xt <- xtabs(~age+chcond, data=ahs)
mosaicplot(xt)

xt <- xtabs(~sex+chcond, data=ahs)
mosaicplot(xt)

##----------------------------------------------------------------------

library(lattice)
library(latticeExtra)

useOuterStrips(
    combineLimits(
        xyplot(Ndoc+Nndoc+Nadm+Nhosp+Nmed~age|sex,
               outer=TRUE, data=ahs,
               jitter.x=TRUE, amount=0.01,
               type=c("p", "a"),
               scales=list(y=list(relation="free")),
               ylab="Number or occurences",
               xlab="Age (years/100)")
    )
)

useOuterStrips(
    combineLimits(
        xyplot(Ndoc+Nndoc+Nadm+Nhosp+Nmed~income|sex,
               outer=TRUE, data=ahs,
               jitter.x=TRUE, amount=0.01,
               type=c("p", "a"),
               scales=list(y=list(relation="free")),
               ylab="Number or occurences",
               xlab="Income")
    )
)

useOuterStrips(
    combineLimits(
        xyplot(Ndoc+Nndoc+Nadm+Nhosp+Nmed~age|chcond,
               groups=sex, outer=TRUE, data=ahs,
               jitter.x=TRUE, amount=0.01,
               type=c("p", "a"),
               scales=list(y=list(relation="free")),
               ylab="Number or occurences",
               xlab="Age (years/100)")
    )
)