diff --git a/.Rbuildignore b/.Rbuildignore index c7b544af4620fe2f6b43e12da6dc0d6ba159b2f3..e3e71d6596a42ae7c51d9f320c68448f61104086 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -8,3 +8,6 @@ buildPkg.R README.Rmd Examples/ +^data-raw$ +\#*\# +^\.\#* diff --git a/.gitignore b/.gitignore index b46f10df4f23a0ce962b548f7e0483182b0c5286..d55718fc065746d8d88dec22a25d5e1f15e1bd52 100644 --- a/.gitignore +++ b/.gitignore @@ -12,4 +12,6 @@ ^\.Rd2pdf* .Rd2pdf5504 .Rd2pdf5516 -.Rd2pdf* +.Rd2pdf* +data-raw/*.txt +inst/doc diff --git a/DESCRIPTION b/DESCRIPTION index 41e758dd1c9718087d5a57ad169db0e7dfd9aac2..b8a49da74eb2bfdc423f63e6bb205cbe44d29bde 100644 --- a/DESCRIPTION +++ b/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 diff --git a/R/mc_S3_methods.R b/R/mc_S3_methods.R index 2219d5f6116c91890b22806bf4461ed691ec0b40..04a6d6b56a2f0b9dd810ba2e09244c17c321ba0e 100644 --- a/R/mc_S3_methods.R +++ b/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,8 +271,9 @@ plot.mcglm <- function(x, type = "residuals", ...) { plot(res ~ fit_values, ylab = "Pearson residuals", xlab = "Fitted values") - temp <- loess.smooth(fitted(object)[, i], - residuals(object, type = "pearson")[, i]) + temp <- loess.smooth( + fitted(object)[, i], + residuals(object, type = "pearson")[, i]) lines(temp$x, temp$y) qqnorm(res) qqline(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") diff --git a/R/mcglm.R b/R/mcglm.R index 84a11369430944b88a905ff6170a00760ed89b36..bd5bb010d760adc3e4fe47ef5f5ea9ed7982227a 100644 --- a/R/mcglm.R +++ b/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 diff --git a/buildPkg.R b/buildPkg.R index 785333db8a80f615dbbfc5d91657206bc1fb3542..e7f60d77c836a39057779fd94e523cd634ca5485 100644 --- a/buildPkg.R +++ b/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. diff --git a/data-raw/RData2txt.R b/data-raw/RData2txt.R new file mode 100644 index 0000000000000000000000000000000000000000..26918341b921bd365394b55db225018ab897b6c8 --- /dev/null +++ b/data-raw/RData2txt.R @@ -0,0 +1,46 @@ +##====================================================================== +## 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) + +##---------------------------------------------------------------------- diff --git a/data-raw/ahs.R b/data-raw/ahs.R new file mode 100644 index 0000000000000000000000000000000000000000..fd48eb6bb15d7d2a7dcadb76395350f65f311aee --- /dev/null +++ b/data-raw/ahs.R @@ -0,0 +1,50 @@ +##---------------------------------------------------------------------- +## 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)") + ) +) + +##---------------------------------------------------------------------- diff --git a/man/ahs.Rd b/man/ahs.Rd index 7cb32800493c4d577077f654dad50863d8004854..7f159c406fc11a19c969027f848e65ad5c5ed678 100644 --- a/man/ahs.Rd +++ b/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} diff --git a/man/plot.mcglm.Rd b/man/plot.mcglm.Rd index ebb64e8fd5d3dc830eeba23dd354d56390fe49d7..d0c270b0fb9c0e03734b6a395d982143d9026123 100644 --- a/man/plot.mcglm.Rd +++ b/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", ...) } diff --git a/vignettes/vignette-01.Rmd b/vignettes/vignette-01.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..9812fb1824aec2ee95588d39fe59b28f16c73796 --- /dev/null +++ b/vignettes/vignette-01.Rmd @@ -0,0 +1,238 @@ +--- +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