Merge branch 'issue#7' into 'devel'

Issue#7: Data set and vignette

This branch:
  - Add `data-raw/` directory. In this directory, data preparation is done. The code to `@examples` section in the documentation is also done there.
  - Add R script to get `txt` files from `RData` files.
  - Add `data/ahs.RData`, a data.frame named `ahs`;
  - Modifies de data.frame and improves de documentation (CMT 0a92ba26 for details);
  - Add lattice and latticeExtra in the Suggests, as they are required in the `@examples` of `ahs`.
  - Add a vignette template that contains plots and summaries of the `ahs` data set.
  - Add build_vignette() in buildPkg.R.
  - Build the package and vignette.
  - Compilation goes without errors. Success!.


See merge request !5

.Rbuildignore

@@ -8,3 +8,6 @@
 buildPkg.R
 README.Rmd
 Examples/
+^data-raw$
+\#*\#
+^\.\#*

.gitignore

@@ -13,3 +13,5 @@
 .Rd2pdf5504
 .Rd2pdf5516
 .Rd2pdf*
+data-raw/*.txt
+inst/doc

DESCRIPTION

@@ -14,7 +14,9 @@ Depends:
     R (>= 3.2.1)
 Suggests:
     testthat,
-    plyr
+    plyr,
+    lattice,
+    latticeExtra
 Imports:
     Matrix,
     assertthat
@@ -23,3 +25,4 @@ LazyData: TRUE
 URL: http://git.leg.ufpr.br/wbonat/mcglm
 BugReports: http://git.leg.ufpr.br/wbonat/mcglm/issues
 Encoding: UTF-8
+VignetteBuilder: knitr

R/mc_S3_methods.R

@@ -239,7 +239,9 @@ fitted.mcglm <- function(object, ...) {
     output <- Matrix(object$fitted, ncol = n_resp, nrow = object$n_obs)
     return(output)
 }
-#' @title Default Multivariate Covariance Generalized Linear models plotting
+
+#' @title Default Multivariate Covariance Generalized Linear models
+#'     plotting
 #' @name plot.mcglm
 #'
 #' @description takes a fitted \code{mcglm} object and do plots based on
@@ -269,7 +271,8 @@ plot.mcglm <- function(x, type = "residuals", ...) {
             plot(res ~ fit_values,
                  ylab = "Pearson residuals",
                  xlab = "Fitted values")
-            temp <- loess.smooth(fitted(object)[, i],
+            temp <- loess.smooth(
+                fitted(object)[, i],
                 residuals(object, type = "pearson")[, i])
             lines(temp$x, temp$y)
             qqnorm(res)
@@ -362,7 +365,8 @@ print.mcglm <- function(x, ...) {
         names(tau_temp) <- rep("", length(tau_temp))
         print(tau_temp)
         cat("\n")
-        power_temp <- coef(object, response = i, type = "power")$Estimate
+        power_temp <- coef(object, response = i,
+                           type = "power")$Estimate
         if (length(power_temp) != 0) {
             names(power_temp) <- ""
             cat("Power:\n")
@@ -445,13 +449,15 @@ summary.mcglm <- function(object, ...) {
         cat("Covariance function:", object$covariance[[i]])
         cat("\n")
         cat("Regression:\n")
-        tab_beta <- coef(object, std.error = TRUE, response = i, type = "beta")[, 1:2]
+        tab_beta <- coef(object, std.error = TRUE,
+                         response = i, type = "beta")[, 1:2]
         tab_beta$"Z value" <- tab_beta[, 1]/tab_beta[, 2]
         rownames(tab_beta) <- object$beta_names[[i]]
         output[i][[1]]$Regression <- tab_beta
         print(tab_beta)
         cat("\n")
-        tab_power <- coef(object, std.error = TRUE, response = i, type = "power")[, 1:2]
+        tab_power <- coef(object, std.error = TRUE,
+                          response = i, type = "power")[, 1:2]
         tab_power$"Z value" <- tab_power[, 1]/tab_power[, 2]
         rownames(tab_power) <- NULL
         if (dim(tab_power)[1] != 0) {
@@ -461,14 +467,16 @@ summary.mcglm <- function(object, ...) {
             cat("\n")
         }
         cat("Dispersion:\n")
-        tab_tau <- coef(object, std.error = TRUE, response = i, type = "tau")[, 1:2]
+        tab_tau <- coef(object, std.error = TRUE,
+                        response = i, type = "tau")[, 1:2]
         tab_tau$"Z value" <- tab_tau[, 1]/tab_tau[, 2]
         rownames(tab_tau) <- NULL
         output[i][[1]]$tau <- tab_tau
         print(tab_tau)
         cat("\n")
     }
-    tab_rho <- coef(object, std.error = TRUE, response = NA, type = "correlation")[, c(3, 1, 2)]
+    tab_rho <- coef(object, std.error = TRUE,
+                    response = NA, type = "correlation")[, c(3, 1, 2)]
     tab_rho$"Z value" <- tab_rho[, 2]/tab_rho[, 3]
     if (dim(tab_rho)[1] != 0) {
         cat("Correlation matrix:\n")

R/mcglm.R

@@ -33,41 +33,40 @@ NULL
 #'
 #' \itemize{
 #'
-#' \item \code{sex} - Factor, two levels (0-Male; 1-Female).
+#' \item \code{sex} - Factor with levels \code{male} and \code{female}.
 #'
 #' \item \code{age} - Respondent's age in years divided by 100.
 #'
 #' \item \code{income} - Respondent's annual income in Australian
 #'     dollars divided by 1000.
 #'
-#' \item \code{levyplus} - Factor, two levels (1- if respondent is
-#'     covered by private health insurance fund for private patients in
-#'     public hospital (with doctor of choice); 0 - otherwise).
+#' \item \code{levyplus} - Coded factor. If respondent is covered by
+#'     private health insurance fund for private patients in public
+#'     hospital with doctor of choice (1) or otherwise (0).
 #'
-#' \item \code{freepoor} - Factor, two levels (1 - if respondent is
-#'     covered by government because low income, recent immigrant,
-#'     unemployed; 0 - otherwise).
+#' \item \code{freepoor} - Coded factor. If respondent is covered by
+#'     government because low income, recent immigrant, unemployed (1)
+#'     or otherwise (0).
 #'
-#' \item \code{freerepa} - Factor, two levels (1 - if respondent is
-#'     covered free by government because of old-age or disability
-#'     pension, or because invalid veteran or family of deceased
-#'     veteran; 0 - otherwise).
+#' \item \code{freerepa} - Coded factor. If respondent is covered free
+#'     by government because of old-age or disability pension, or
+#'     because invalid veteran or family of deceased veteran (1) or
+#'     otherwise (0).
 #'
 #' \item \code{illnes} - Number of illnesses in past 2 weeks, with 5 or
-#'     more weeks coded as 5.
+#'     illnesses coded as 5.
 #'
 #' \item \code{actdays} - Number of days of reduced activity in the past
 #'     two weeks due to illness or injury.
 #'
 #' \item \code{hscore} - Respondent's general health questionnaire score
-#'     using Goldberg's method; high score indicates poor health.
+#'     using Goldberg's method. High score indicates poor health.
 #'
-#' \item \code{chcond1} - Factor, two levels (1 - if respondent has
-#'     chronic condition(s) but is not limited in activity; 0 -
-#'     otherwise).
-#'
-#' \item \code{chcond2} - Factor, two levels (1 if respondent has
-#'     chronic condition(s) and is limited in activity; 0 - otherwise).
+#' \item \code{chcond} - Factor with three levels. If respondent has
+#'     chronic condition(s) and is limited in activity (\code{limited}),
+#'     or if the respondent has chronic condition(s) but is not limited
+#'     in activity (\code{nonlimited}) or otherwise (\code{otherwise},
+#'     reference level).
 #'
 #' \item \code{Ndoc} - Number of consultations with a doctor or
 #'     specialist (response variable).
@@ -85,8 +84,6 @@ NULL
 #' \item \code{Nmed} - Total number of prescribed and non prescribed
 #'     medications used in the past two days.
 #'
-#' \item \code{id} - Respondent's index.
-#'
 #' }
 #'
 #' @docType data
@@ -101,4 +98,42 @@ NULL
 #'     the elderly: A finite mixture approach, Journal of Applied
 #'     Econometrics 12(3):313--336.
 #'
+#' @examples
+#'
+#' library(lattice)
+#' library(latticeExtra)
+#'
+#' data(ahs, package="mcglm")
+#' str(ahs)
+#'
+#' xt <- xtabs(~age+sex, data=ahs)
+#' mosaicplot(xt)
+#'
+#' xt <- xtabs(~age+chcond, data=ahs)
+#' mosaicplot(xt)
+#'
+#' 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="Age (years/100)")
+#'     )
+#' )
+#'
 NULL

buildPkg.R

@@ -2,7 +2,11 @@
 ## Script to build and verify the package.
 
 if(!grepl(x=getwd(), pattern="/mcglm$")){
-    stop("Move to /mcglm directory.")
+    if (Sys.info()["user"]=="walmes"){
+        setwd("~/GitLab/mcglm")
+    }
+    ## stop("Move to /mcglm directory.")
+    cat(getwd(), "\n")
 }
 
 ##----------------------------------------------------------------------
@@ -12,7 +16,6 @@ library(devtools)
 
 ## Load the package (to make functiona available).
 load_all()
-search()
 
 ## Create/update NAMESPACE, *.Rd files.
 document()
@@ -39,6 +42,18 @@ build(manual = TRUE, vignettes = FALSE)
 # build the binary version for windows (not used)
 # build_win()
 
+##----------------------------------------------------------------------
+## Package vignette.
+## Based on: http://r-pkgs.had.co.nz/vignettes.html
+
+## Create the vignette template. Do just once.
+## use_vignette("vignette-01")
+
+build_vignettes()
+
+## vignette()
+## vignette("vignette-01", package="mcglm")
+
 ##----------------------------------------------------------------------
 ## Generate the README.md.
 

data-raw/RData2txt.R

+##======================================================================
+## Generate plain text files from the RData files for all datasets in
+## the package.
+##======================================================================
+
+getwd()
+f <- list.files(path="../data", pattern="*.RData")
+
+sapply(f,
+       FUN=function(f){
+           load(paste0("../data/", f))
+           g <- sub("\\.RData", "", f)
+           txt <- paste0(g, ".txt")
+           dataset <- eval(parse(text=g))
+           cat(file=txt, sep="\n",
+               "##----------------------------------------------------------------------",
+               "## This dataset is part of mcglm package.",
+               "## Visit http://git.leg.ufpr.br/wbonat/mcglm for details.",
+               "##----------------------------------------------------------------------")
+           write.table(x=dataset,
+                       file=txt,
+                       sep="\t", quote=FALSE, row.names=FALSE,
+                       append=TRUE)
+       })
+
+##----------------------------------------------------------------------
+## Uploading files to the public folder: www.leg.ufpr.br/~leg/mcglm.
+
+## Port and IP.
+u <- scan(n=2, what=character())
+
+## Verifica o conteúdo do diretório /datasets.
+cmd <- paste0("ssh leg@", u[2], " -p", u[1],
+              " \"ls -la ~/public_html/mcglm/datasets\"")
+system(cmd)
+
+## By rsync.
+cmd <- paste0("rsync -avzh --progress *.txt ",
+              "-e \"ssh -p ", u[1], "\" leg@", u[2],
+              ":~/public_html/mcglm/datasets/")
+system(cmd)
+
+url <- "http://www.leg.ufpr.br/~leg/mcglm/datasets"
+browseURL(url)
+
+##----------------------------------------------------------------------

data-raw/ahs.R

+##----------------------------------------------------------------------
+## Prepare de the data set.
+
+setwd("/home/walmes/GitLab/mcglm/data-raw")
+
+ahs <- read.table("ahs.txt", header=TRUE, sep="\t")
+str(ahs)
+
+## save(ahs, file="../data/ahs.RData")
+
+##----------------------------------------------------------------------
+## Include in the @examples.
+
+library(lattice)
+library(latticeExtra)
+
+data(ahs, package="mcglm")
+str(ahs)
+
+xt <- xtabs(~age+sex, data=ahs)
+mosaicplot(xt)
+
+xt <- xtabs(~age+chcond, data=ahs)
+mosaicplot(xt)
+
+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="Age (years/100)")
+    )
+)
+
+##----------------------------------------------------------------------

man/ahs.Rd

@@ -22,41 +22,40 @@ The Australian health survey was used by Bonat and
 
 \itemize{
 
-\item \code{sex} - Factor, two levels (0-Male; 1-Female).
+\item \code{sex} - Factor with levels \code{male} and \code{female}.
 
 \item \code{age} - Respondent's age in years divided by 100.
 
 \item \code{income} - Respondent's annual income in Australian
     dollars divided by 1000.
 
-\item \code{levyplus} - Factor, two levels (1- if respondent is
-    covered by private health insurance fund for private patients in
-    public hospital (with doctor of choice); 0 - otherwise).
+\item \code{levyplus} - Coded factor. If respondent is covered by
+    private health insurance fund for private patients in public
+    hospital with doctor of choice (1) or otherwise (0).
 
-\item \code{freepoor} - Factor, two levels (1 - if respondent is
-    covered by government because low income, recent immigrant,
-    unemployed; 0 - otherwise).
+\item \code{freepoor} - Coded factor. If respondent is covered by
+    government because low income, recent immigrant, unemployed (1)
+    or otherwise (0).
 
-\item \code{freerepa} - Factor, two levels (1 - if respondent is
-    covered free by government because of old-age or disability
-    pension, or because invalid veteran or family of deceased
-    veteran; 0 - otherwise).
+\item \code{freerepa} - Coded factor. If respondent is covered free
+    by government because of old-age or disability pension, or
+    because invalid veteran or family of deceased veteran (1) or
+    otherwise (0).
 
 \item \code{illnes} - Number of illnesses in past 2 weeks, with 5 or
-    more weeks coded as 5.
+    illnesses coded as 5.
 
 \item \code{actdays} - Number of days of reduced activity in the past
     two weeks due to illness or injury.
 
 \item \code{hscore} - Respondent's general health questionnaire score
-    using Goldberg's method; high score indicates poor health.
+    using Goldberg's method. High score indicates poor health.
 
-\item \code{chcond1} - Factor, two levels (1 - if respondent has
-    chronic condition(s) but is not limited in activity; 0 -
-    otherwise).
-
-\item \code{chcond2} - Factor, two levels (1 if respondent has
-    chronic condition(s) and is limited in activity; 0 - otherwise).
+\item \code{chcond} - Factor with three levels. If respondent has
+    chronic condition(s) and is limited in activity (\code{limited}),
+    or if the respondent has chronic condition(s) but is not limited
+    in activity (\code{nonlimited}) or otherwise (\code{otherwise},
+    reference level).
 
 \item \code{Ndoc} - Number of consultations with a doctor or
     specialist (response variable).
@@ -74,9 +73,44 @@ The Australian health survey was used by Bonat and
 \item \code{Nmed} - Total number of prescribed and non prescribed
     medications used in the past two days.
 
-\item \code{id} - Respondent's index.
-
 }
 }
+\examples{
+library(lattice)
+library(latticeExtra)
+
+data(ahs, package="mcglm")
+str(ahs)
+
+xt <- xtabs(~age+sex, data=ahs)
+mosaicplot(xt)
+
+xt <- xtabs(~age+chcond, data=ahs)
+mosaicplot(xt)
+
+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="Age (years/100)")
+    )
+)
+}
 \keyword{datasets}
 

man/plot.mcglm.Rd

@@ -2,7 +2,8 @@
 % Please edit documentation in R/mc_S3_methods.R
 \name{plot.mcglm}
 \alias{plot.mcglm}
-\title{Default Multivariate Covariance Generalized Linear models plotting}
+\title{Default Multivariate Covariance Generalized Linear models
+    plotting}
 \usage{
 \method{plot}{mcglm}(x, type = "residuals", ...)
 }

vignettes/vignette-01.Rmd

+---
+title: "Multivariate Covariance Generalized Linear Models"
+author: "Wagner Hugo Bonat"
+date: "`r paste('mcglm', packageVersion('mcglm'), Sys.Date())`"
+output:
+  rmarkdown::html_vignette:
+    fig_width: 6
+    fig_height: 6
+    toc: true
+    toc_dep: 3
+vignette: >
+  %\VignetteIndexEntry{Multivariate Covariance Generalized Linear Models}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+<style type="text/css">
+body, td, caption {
+    font-family: "Palatino Linotype", "Book Antiqua", Palatino, serif;
+    background-color: white;
+    font-size: 16px;
+}
+
+tt, code, pre {
+    font-family: "Inconsolata", "Andale Mono", monospace;
+}
+
+code {
+    font-size: 16px;
+}
+
+pre code {
+    font-size: 14px;
+}
+
+pre:not([class]) code {
+    background-color: #92BFB1;
+}
+pre, code {
+    background-color: #62BFB1;
+    border-radius: 3px;
+    color: #333;
+}
+
+/* R output */
+pre:not([class]) code {
+    background-color: #D4D4D4;
+}
+pre:not([class]), code {
+    background-color: #D4D4D4;
+}
+
+/* R input */
+pre, code {
+    border-radius: 3px;
+    background-color: #EDEDED;
+    color: #333;
+}
+
+img {
+    max-width: 100% !important;
+    display: block;
+    margin: auto;
+}
+
+.MathJax {
+    font-size: 80% !important;
+}
+
+</style>
+
+```{r setup, include=FALSE}
+##----------------------------------------------------------------------
+
+library(knitr)
+
+opts_chunk$set(
+    dev.args=list(family="Palatino"))
+
+options(width=68)
+
+##----------------------------------------------------------------------
+
+library(latticeExtra)
+rm(list=ls())
+
+## Color palette.
+mycol <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00",
+           "#FFFF33")
+dput(mycol)
+
+## Trellis graphical style.
+ps <- list(
+    box.rectangle=list(col=1, fill=c("gray70")),
+    box.umbrella=list(col=1, lty=1),
+    dot.symbol=list(col=1, pch=19),
+    dot.line=list(col="gray50", lty=3),
+    plot.symbol=list(col=1, cex=0.8),
+    plot.line=list(col=1),
+    plot.polygon=list(col="gray95"),
+    superpose.line=list(col=mycol, lty=1),
+    superpose.symbol=list(col=mycol, pch=1),
+    superpose.polygon=list(col=mycol),
+    strip.background=list(col=c("gray80","gray50"))
+    )
+trellis.par.set(ps)
+## show.settings()
+
+```
+
+****
+
+To install the stable version of [`mcglm`][], use
+`devtools::install_git()`. For more information, visit [mcglm/README].
+
+```{r, eval=FALSE}
+library(devtools)
+install_git("http://git.leg.ufpr.br/wbonat/mcglm.git")
+```
+
+```{r, eval=FALSE, error=FALSE, message=FALSE, warning=FALSE}
+library(mcglm)
+packageVersion("mcglm")
+```
+
+```{r, echo=FALSE, error=FALSE, message=FALSE, warning=FALSE}
+library(mcglm)
+packageVersion("mcglm")
+```
+
+****
+## The Australian Health Survey
+
+```{r}
+##----------------------------------------------------------------------
+## Loadin the Australian Health Survey data.
+
+data(ahs)
+
+## Object structure.
+str(ahs)
+
+## Descriptive measures.
+summary(ahs)
+
+##----------------------------------------------------------------------
+## Frequency tables.
+
+names(ahs)[c(1, 4:7, 10)]
+
+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)")
+    )
+)
+
+```
+
+<!---------------------------------------------------------------------- -->
+
+[`mcglm`]: http://git.leg.ufpr.br/wbonat/mcglm
+[mcglm/README]: http://git.leg.ufpr.br/wbonat/mcglm/blob/master/README.md