merge conflict

NEWS.md

@@ -2,6 +2,12 @@
 
 ## merTools 1.1
 
+### Bug fixes
+
+### New Functionality
+
+- Add support for `probit` models and limited support for other `glmm` link functions, with warning
+- Add ability for user-specified seed for reproducibility
 - Add support for `blmer` objects from the `blme` package
 
 ## merTools 0.1

R/merData.R

@@ -39,12 +39,13 @@ stripAttributes <- function(data){
 #' @param type what kind of draw to make. Options include random or average
 #' @param varList a list specifying filters to subset the data by when making the
 #' draw
+#' @param seed numeric, optional argument to set seed for simulations, ignored if type="average"
 #' @return a data.frame with a single row representing the desired observation
 #' @details In cases of tie, ".", may be substituted for factors.
 #' @export draw
 #' @rdname draw
 draw <- function(object, type = c("random", "average"),
-                 varList = NULL){
+                 varList = NULL, seed = NULL){
   UseMethod("draw")
 }
 
@@ -63,10 +64,10 @@ draw <- function(object, type = c("random", "average"),
 #' draw(fm1, type = "average", varList = list("Subject" = "308"))
 #'
 draw.merMod <- function(object, type = c("random", "average"),
-                        varList = NULL){
+                        varList = NULL, seed = NULL){
   type <- match.arg(type, c("random", "average"), several.ok = FALSE)
   if(type == 'random'){
-    out <- randomObs(object, varList)
+    out <- randomObs(object, varList, seed)
     return(out)
   } else if(type == 'average'){
     out <- averageObs(object, varList)
@@ -79,14 +80,21 @@ draw.merMod <- function(object, type = c("random", "average"),
 #' @description Select a random observation from the model frame of a merMod
 #' @param merMod an object of class merMod
 #' @param varList optional, a named list of conditions to subset the data on
+#' @param seed numeric, optional argument to set seed for simulations
 #' @return a data frame with a single row for a random observation, but with full
 #' factor levels. See details for more.
 #' @details Each factor variable in the data frame has all factor levels from the
 #' full model.frame stored so that the new data is compatible with predict.merMod
-randomObs <- function(merMod, varList){
+randomObs <- function(merMod, varList, seed = NULL){
   if(!missing(varList)){
     data <- subsetList(merMod@frame, varList)
   }
+
+  if (!is.null(seed))
+    set.seed(seed)
+  else if (!exists(".Random.seed", envir = .GlobalEnv))
+    runif(1)
+
   out <- data[sample(1:nrow(data), 1),]
   chars <- !sapply(out, is.numeric)
   for(i in names(out[, chars])){

R/merExtract.R

@@ -53,6 +53,7 @@ REextract <- function(merMod){
 #' @param merMod a merMod object from the lme4 package
 #' @param n.sims number of simulations to use
 #' @param oddsRatio logical, should parameters be converted to odds ratios?
+#' @param seed numeric, optional argument to set seed for simulations
 #' @importFrom arm sim
 #' @import lme4
 #' @return a data frame with the following columns
@@ -72,9 +73,14 @@ REextract <- function(merMod){
 #' re2 <- REsim(m2, 25)
 #' head(re2)
 #' @export
-REsim <- function(merMod, n.sims = 200, oddsRatio = FALSE){
+REsim <- function(merMod, n.sims = 200, oddsRatio = FALSE, seed=NULL){
   stopifnot(class(merMod) %in% c("lmerMod", "glmerMod", "blmerMod",
                                  "bglmerMod"))
+  if (!is.null(seed))
+    set.seed(seed)
+  else if (!exists(".Random.seed", envir = .GlobalEnv))
+    runif(1)
+
   mysim <- arm::sim(merMod, n.sims = n.sims)
   reDims <- length(mysim@ranef)
   tmp.out <- vector("list", reDims)
@@ -105,6 +111,7 @@ REsim <- function(merMod, n.sims = 200, oddsRatio = FALSE){
 #' @param merMod a merMod object from the lme4 package
 #' @param n.sims number of simulations to use
 #' @param oddsRatio logical, should parameters be converted to odds ratios?
+#' @param seed numeric, optional argument to set seed for simulations
 #' @importFrom arm sim
 #' @import lme4
 #' @return a data frame with the following columns
@@ -122,9 +129,14 @@ REsim <- function(merMod, n.sims = 200, oddsRatio = FALSE){
 #' fe2 <- FEsim(m2, 25)
 #' head(fe2)
 #' @export
-FEsim <- function(merMod, n.sims = 200, oddsRatio=FALSE){
+FEsim <- function(merMod, n.sims = 200, oddsRatio=FALSE, seed=NULL){
   stopifnot(class(merMod) %in% c("lmerMod", "glmerMod", "blmerMod",
                                  "bglmerMod"))
+  if (!is.null(seed))
+    set.seed(seed)
+  else if (!exists(".Random.seed", envir = .GlobalEnv))
+    runif(1)
+
   mysim <- arm::sim(merMod, n.sims = n.sims)
   means <- apply(mysim@fixef, MARGIN = 2, mean)
   medians <- apply(mysim@fixef, MARGIN = 2, median)

R/merPredict.R

@@ -15,6 +15,7 @@
 #' @param include.resid.var logical, include or exclude the residual variance for
 #' linear models
 #' @param returnSims logical, should all n.sims simulations be returned?
+#' @param seed numeric, optional argument to set seed for simulations
 #' @return a data.frame iwth three columns:
 #' \describe{
 #'     \item{\code{fit}}{The center of the distribution of predicted values as defined by
@@ -71,7 +72,8 @@
 predictInterval <- function(merMod, newdata, level = 0.95,
                             n.sims=100, stat=c("median","mean"),
                             type=c("linear.prediction", "probability"),
-                            include.resid.var=TRUE, returnSims = FALSE){
+                            include.resid.var=TRUE, returnSims = FALSE,
+                            seed=NULL){
   outs <- newdata
   predict.type <- match.arg(type,
                             c("linear.prediction", "probability"),
@@ -80,6 +82,11 @@ predictInterval <- function(merMod, newdata, level = 0.95,
                          c("median","mean"),
                          several.ok = FALSE)
 
+  if (!is.null(seed))
+    set.seed(seed)
+  else if (!exists(".Random.seed", envir = .GlobalEnv))
+    runif(1)
+
   ##First: check if it is a GLMM or NLMM and draw from sigma distribution or incorporate scale parameter if GLMM
   merMod.devcomp <- getME(merMod, "devcomp")
   if (merMod.devcomp$dims[["GLMM"]] == 0 &

R/subBoot.R

@@ -22,7 +22,7 @@ thetaExtract <- function(merMod){
 #' in the merMod object, by default will resample the entire model frame
 #' @param FUN the function to apply to each bootstrapped model
 #' @param R the number of bootstrap replicates, default is 100
-#'
+#' @param seed numeric, optional argument to set seed for simulations
 #' @return a data.frame of parameters extracted from each of the R replications.
 #' The original values are appended to the top of the matrix.
 #' @details This function allows users to estimate parameters of a
@@ -33,12 +33,18 @@ thetaExtract <- function(merMod){
 #' resultMatrix <- subBoot(fm1, n = 160, FUN = thetaExtract, R = 20)
 #' }
 #' @export
-subBoot <- function(merMod, n = NULL, FUN, R = 100){
+subBoot <- function(merMod, n = NULL, FUN, R = 100, seed=NULL){
   if(missing(n))
     n <- nrow(merMod@frame)
   resultMat <- matrix(FUN(merMod), nrow = 1)
   tmp <- matrix(data=NA, nrow=R, ncol=ncol(resultMat))
   resultMat <- rbind(resultMat, tmp); rm(tmp)
+
+  if (!is.null(seed))
+    set.seed(seed)
+  else if (!exists(".Random.seed", envir = .GlobalEnv))
+    runif(1)
+
     for(i in 1:R){
       rows <- as.numeric(row.names(merMod@frame))
       mysamp <- as.character(sample(rows,  n, replace=TRUE))

README.Rmd

@@ -24,7 +24,7 @@ by Jared E. Knowles and Carl Frederick
 [![Coverage Status](https://coveralls.io/repos/jknowles/merTools/badge.svg?branch=master)](https://coveralls.io/r/jknowles/merTools?branch=master)
 [![Github Issues](http://githubbadges.herokuapp.com/jknowles/merTools/issues.svg)](https://github.com/jknowles/merTools/issues)
 [![Pending Pull-Requests](http://githubbadges.herokuapp.com/jknowles/merTools/pulls.svg?style=flat)](https://github.com/jknowles/merTools/pulls)
-[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/caretEnsemble)](http://cran.r-project.org/web/packages/merTools)
+[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/merTools)](http://cran.r-project.org/web/packages/merTools)
 [![Downloads](http://cranlogs.r-pkg.org/badges/merTools)](http://cran.rstudio.com/package=merTools)
 
 Working with generalized linear mixed models (GLMM) and linear mixed models (LMM) 
@@ -43,7 +43,7 @@ those tools.
 library(devtools)
 install_github("jknowles/merTools")
 
-# CRAN version -- coming soon
+# CRAN version
 install.packages("merTools")
 ```
 

README.md

@@ -6,7 +6,7 @@ A package for getting the most of our multilevel models in R
 
 by Jared E. Knowles and Carl Frederick
 
-[![Travis-CI Build Status](https://travis-ci.org/jknowles/merTools.png?branch=master)](https://travis-ci.org/jknowles/merTools) [![Coverage Status](https://coveralls.io/repos/jknowles/merTools/badge.svg?branch=master)](https://coveralls.io/r/jknowles/merTools?branch=master) [![Github Issues](http://githubbadges.herokuapp.com/jknowles/merTools/issues.svg)](https://github.com/jknowles/merTools/issues) [![Pending Pull-Requests](http://githubbadges.herokuapp.com/jknowles/merTools/pulls.svg?style=flat)](https://github.com/jknowles/merTools/pulls) [![CRAN\_Status\_Badge](http://www.r-pkg.org/badges/version/caretEnsemble)](http://cran.r-project.org/web/packages/merTools) [![Downloads](http://cranlogs.r-pkg.org/badges/merTools)](http://cran.rstudio.com/package=merTools)
+[![Travis-CI Build Status](https://travis-ci.org/jknowles/merTools.png?branch=master)](https://travis-ci.org/jknowles/merTools) [![Coverage Status](https://coveralls.io/repos/jknowles/merTools/badge.svg?branch=master)](https://coveralls.io/r/jknowles/merTools?branch=master) [![Github Issues](http://githubbadges.herokuapp.com/jknowles/merTools/issues.svg)](https://github.com/jknowles/merTools/issues) [![Pending Pull-Requests](http://githubbadges.herokuapp.com/jknowles/merTools/pulls.svg?style=flat)](https://github.com/jknowles/merTools/pulls) [![CRAN\_Status\_Badge](http://www.r-pkg.org/badges/version/merTools)](http://cran.r-project.org/web/packages/merTools) [![Downloads](http://cranlogs.r-pkg.org/badges/merTools)](http://cran.rstudio.com/package=merTools)
 
 Working with generalized linear mixed models (GLMM) and linear mixed models (LMM) has become increasingly easy with advances in the `lme4` package. As we have found ourselves using these models more and more within our work, we, the authors, have developed a set of tools for simplifying and speeding up common tasks for interacting with `merMod` objects from `lme4`. This package provides those tools.
 
@@ -18,7 +18,7 @@ Installation
 library(devtools)
 install_github("jknowles/merTools")
 
-# CRAN version -- coming soon
+# CRAN version
 install.packages("merTools")
 ```
 
@@ -66,16 +66,16 @@ With `predictInterval` we obtain predictions that are more like the standard obj
 predictInterval(m1, newdata = InstEval[1:10, ], n.sims = 500, level = 0.9, 
                 stat = 'median')
 #>         fit      lwr      upr
-#> 1  3.119054 1.065297 5.269603
-#> 2  3.188968 1.151337 4.996861
-#> 3  3.355738 1.412177 5.477573
-#> 4  3.072446 1.069449 5.204614
-#> 5  3.247946 1.294070 5.280065
-#> 6  3.255983 1.330520 5.095135
-#> 7  4.184949 2.311159 6.095691
-#> 8  3.803738 1.912956 5.720153
-#> 9  3.707334 1.621160 5.816649
-#> 10 3.415814 1.530551 5.388098
+#> 1  3.205852 1.119040 5.073677
+#> 2  3.014957 1.188017 4.894843
+#> 3  3.540154 1.853743 5.582607
+#> 4  3.068565 1.140955 5.257027
+#> 5  3.106084 1.423745 5.191017
+#> 6  3.112947 1.287469 5.414535
+#> 7  4.204398 2.183662 6.274698
+#> 8  3.894898 1.784715 5.885646
+#> 9  3.795654 1.757574 6.010877
+#> 10 3.311081 1.307305 5.291231
 ```
 
 Note that `predictInterval` is slower because it is computing simulations. It can also return all of the simulated `yhat` values as an attribute to the predict object itself.
@@ -91,12 +91,12 @@ Plotting
 feSims <- FEsim(m1, n.sims = 100)
 head(feSims)
 #>          term        mean      median         sd
-#> 1 (Intercept)  3.22591531  3.22733117 0.02090560
-#> 2    service1 -0.07157777 -0.07205579 0.01359473
-#> 3   lectage.L -0.18571387 -0.18647789 0.01571551
-#> 4   lectage.Q  0.02497581  0.02385362 0.01305723
-#> 5   lectage.C -0.02565649 -0.02649764 0.01459933
-#> 6   lectage^4 -0.02122105 -0.02163194 0.01360913
+#> 1 (Intercept)  3.22487833  3.22311114 0.01712966
+#> 2    service1 -0.06723781 -0.06652520 0.01269254
+#> 3   lectage.L -0.18628490 -0.18722651 0.01534445
+#> 4   lectage.Q  0.02358642  0.02403768 0.01270285
+#> 5   lectage.C -0.02486989 -0.02294295 0.01314396
+#> 6   lectage^4 -0.01734253 -0.01696322 0.01410475
 ```
 
 And we can also plot this:
@@ -113,12 +113,12 @@ We can also quickly make caterpillar plots for the random-effect terms:
 reSims <- REsim(m1, n.sims = 100)
 head(reSims)
 #>   groupFctr groupID        term        mean      median        sd
-#> 1         s       1 (Intercept)  0.11406193  0.15386504 0.3086416
-#> 2         s       2 (Intercept) -0.05194233 -0.07755123 0.3451832
-#> 3         s       3 (Intercept)  0.27888402  0.27812495 0.3214322
-#> 4         s       4 (Intercept)  0.26677288  0.28692168 0.2801052
-#> 5         s       5 (Intercept)  0.04262341  0.04945323 0.3543950
-#> 6         s       6 (Intercept)  0.11616528  0.10719162 0.2488089
+#> 1         s       1 (Intercept)  0.13651949  0.12518917 0.3028107
+#> 2         s       2 (Intercept) -0.06930294 -0.07109045 0.3048142
+#> 3         s       3 (Intercept)  0.30419535  0.29992844 0.2935636
+#> 4         s       4 (Intercept)  0.33169961  0.34054570 0.2880752
+#> 5         s       5 (Intercept)  0.07810358  0.06039528 0.3427175
+#> 6         s       6 (Intercept)  0.12541267  0.12544140 0.2388220
 ```
 
 ``` r
@@ -153,11 +153,11 @@ impSim <- REimpact(m1, InstEval[7, ], groupFctr = "d", breaks = 5,
                    n.sims = 300, level = 0.9)
 impSim
 #>   case bin   AvgFit     AvgFitSE nobs
-#> 1    1   1 2.799210 2.829958e-04  193
-#> 2    1   2 3.271602 6.527264e-05  240
-#> 3    1   3 3.561880 5.406849e-05  254
-#> 4    1   4 3.864439 5.966251e-05  265
-#> 5    1   5 4.231175 2.176761e-04  176
+#> 1    1   1 2.803774 2.931571e-04  193
+#> 2    1   2 3.262422 6.069156e-05  240
+#> 3    1   3 3.558152 5.751638e-05  254
+#> 4    1   4 3.840147 5.743902e-05  265
+#> 5    1   5 4.239925 1.959778e-04  176
 ```
 
 The result of `REimpact` shows the change in the `yhat` as the case we supplied to `newdata` is moved from the first to the fifth quintile in terms of the magnitude of the group factor coefficient. We can see here that the individual professor effect has a strong impact on the outcome variable. This can be shown graphically as well:

man/FEsim.Rd

@@ -5,14 +5,16 @@
 \title{Simulate fixed effects from merMod
 \code{FEsim} simulates fixed effects from merMod object posterior distributions}
 \usage{
-FEsim(merMod, n.sims = 200, oddsRatio = FALSE)
+FEsim(merMod, n.sims = 200, oddsRatio = FALSE, seed = NULL)
 }
 \arguments{
 \item{merMod}{a merMod object from the lme4 package}
 
 \item{n.sims}{number of simulations to use}
 
 \item{oddsRatio}{logical, should parameters be converted to odds ratios?}
+
+\item{seed}{numeric, optional argument to set seed for simulations}
 }
 \value{
 a data frame with the following columns

man/REsim.Rd

@@ -5,14 +5,16 @@
 \title{Simulate random effects from merMod
 \code{REsim} simulates random effects from merMod object posterior distributions}
 \usage{
-REsim(merMod, n.sims = 200, oddsRatio = FALSE)
+REsim(merMod, n.sims = 200, oddsRatio = FALSE, seed = NULL)
 }
 \arguments{
 \item{merMod}{a merMod object from the lme4 package}
 
 \item{n.sims}{number of simulations to use}
 
 \item{oddsRatio}{logical, should parameters be converted to odds ratios?}
+
+\item{seed}{numeric, optional argument to set seed for simulations}
 }
 \value{
 a data frame with the following columns

man/draw.Rd

@@ -5,9 +5,10 @@
 \alias{draw.merMod}
 \title{Draw a single observation out of an object matching some criteria}
 \usage{
-draw(object, type = c("random", "average"), varList = NULL)
+draw(object, type = c("random", "average"), varList = NULL, seed = NULL)
 
-\method{draw}{merMod}(object, type = c("random", "average"), varList = NULL)
+\method{draw}{merMod}(object, type = c("random", "average"), varList = NULL,
+  seed = NULL)
 }
 \arguments{
 \item{object}{the object to draw from}
@@ -16,6 +17,8 @@ draw(object, type = c("random", "average"), varList = NULL)
 
 \item{varList}{a list specifying filters to subset the data by when making the
 draw}
+
+\item{seed}{numeric, optional argument to set seed for simulations, ignored if type="average"}
 }
 \value{
 a data.frame with a single row representing the desired observation

man/predictInterval.Rd

@@ -6,7 +6,7 @@
 \usage{
 predictInterval(merMod, newdata, level = 0.95, n.sims = 100,
   stat = c("median", "mean"), type = c("linear.prediction", "probability"),
-  include.resid.var = TRUE, returnSims = FALSE)
+  include.resid.var = TRUE, returnSims = FALSE, seed = NULL)
 }
 \arguments{
 \item{merMod}{a merMod object from lme4}
@@ -25,6 +25,8 @@ predictInterval(merMod, newdata, level = 0.95, n.sims = 100,
 linear models}
 
 \item{returnSims}{logical, should all n.sims simulations be returned?}
+
+\item{seed}{numeric, optional argument to set seed for simulations}
 }
 \value{
 a data.frame iwth three columns:

man/randomObs.Rd

@@ -4,12 +4,14 @@
 \alias{randomObs}
 \title{Select a random observation from model data}
 \usage{
-randomObs(merMod, varList)
+randomObs(merMod, varList, seed = NULL)
 }
 \arguments{
 \item{merMod}{an object of class merMod}
 
 \item{varList}{optional, a named list of conditions to subset the data on}
+
+\item{seed}{numeric, optional argument to set seed for simulations}
 }
 \value{
 a data frame with a single row for a random observation, but with full

man/subBoot.Rd

@@ -4,7 +4,7 @@
 \alias{subBoot}
 \title{Bootstrap a subset of an lme4 model}
 \usage{
-subBoot(merMod, n = NULL, FUN, R = 100)
+subBoot(merMod, n = NULL, FUN, R = 100, seed = NULL)
 }
 \arguments{
 \item{merMod}{a valid merMod object}
@@ -15,6 +15,8 @@ in the merMod object, by default will resample the entire model frame}
 \item{FUN}{the function to apply to each bootstrapped model}
 
 \item{R}{the number of bootstrap replicates, default is 100}
+
+\item{seed}{numeric, optional argument to set seed for simulations}
 }
 \value{
 a data.frame of parameters extracted from each of the R replications.