The function `vpc` is failing to stratify on 2 variables.

[Eliford]

I was re-doing a VPC using the vpc package. It used to work when stratifying by 2 variables, but now it doesnot. Below is an rmarkdown file that I use the example given in the vpc package documentation. It works When I use 1 variable to stratify but doesn't with 2 variables. The error message is Error in grep("row", facet) : object 'facet' not found

```{r}
## Load the theophylline PK dataset
library(tidyverse)
library(vpc)
obs <- Theoph
colnames(obs) <- c("id", "wt", "dose", "time", "dv")
obs <- obs %>%
        group_by(id) %>%
        mutate(sex = round(runif(1))) # randomly assign a "sex" covariate
sim <- sim_data(obs, # the design of the dataset
                model = function(x) { # the model
                  pk_oral_1cmt (t = x$time, dose=x$dose * x$wt, ka = x$ka, ke = x$ke, cl = x$cl * x$wt)
                },
                error = list(additive = 0.1),
                theta = c(2.774, 0.0718, .0361),                 # parameter values
                omega_mat = c(0.08854,                           # specified as lower triangle by default;
                              0.02421, 0.02241,                  # note: assumed that every theta has iiv, set to 0 if no iiv.
                              0.008069, 0.008639, 0.02862),
                par_names = c("ka", "ke", "cl"),                 # link the parameters in the model to the thetas/omegas
                n = 500)

#Stratify with 1 variable

obs<-obs%>%mutate(doseg=ifelse(dose<4.5,1,2))
sim<-sim%>%mutate(doseg=ifelse(dose<4.5,1,2))

vpc(sim = sim,
       obs = obs,                               # supply simulation and observation dataframes
       obs_cols = list(
         dv = "dv",                             # these column names are the default,
         idv = "time"),                         # update these if different.
       sim_cols = list(
         dv = "sdv",
         idv = "time"),
       bins = c(0, 2, 4, 6, 8, 10, 16, 25),     # specify bin separators manually
       stratify =c("doseg"),                     # multiple stratifications possible, just supply as vector
       pi = c(0.05, 0.95),                      # prediction interval simulated data to show
       ci = c(0.05, 0.95),                      # confidence intervals to show
       pred_corr = FALSE,                       # perform prediction-correction?
       show = list(obs_dv = TRUE),              # plot observations?
       facet = "wrap",                          # wrap stratifications, or as "row" or "column"
       ylab = "Concentration",
       xlab = "Time (hrs)")

vpc(sim = sim,
       obs = obs,                               # supply simulation and observation dataframes
       obs_cols = list(
         dv = "dv",                             # these column names are the default,
         idv = "time"),                         # update these if different.
       sim_cols = list(
         dv = "sdv",
         idv = "time"),
       bins = c(0, 2, 4, 6, 8, 10, 16, 25),     # specify bin separators manually
       stratify =c("sex"),                     # multiple stratifications possible, just supply as vector
       pi = c(0.05, 0.95),                      # prediction interval simulated data to show
       ci = c(0.05, 0.95),                      # confidence intervals to show
       pred_corr = FALSE,                       # perform prediction-correction?
       show = list(obs_dv = TRUE),              # plot observations?
       facet = "wrap",                          # wrap stratifications, or as "row" or "column"
       ylab = "Concentration",
       xlab = "Time (hrs)")

#Stratify with 2 variables

vpc(sim = sim,
       obs = obs,                               # supply simulation and observation dataframes
       obs_cols = list(
         dv = "dv",                             # these column names are the default,
         idv = "time"),                         # update these if different.
       sim_cols = list(
         dv = "sdv",
         idv = "time"),
       bins = c(0, 2, 4, 6, 8, 10, 16, 25),     # specify bin separators manually
       stratify =c("sex", "doseg"),                     # multiple stratifications possible, just supply as vector
       pi = c(0.05, 0.95),                      # prediction interval simulated data to show
       ci = c(0.05, 0.95),                      # confidence intervals to show
       pred_corr = FALSE,                       # perform prediction-correction?
       show = list(obs_dv = TRUE),              # plot observations?
       facet = "wrap",                          # wrap stratifications, or as "row" or "column"
       ylab = "Concentration",
       xlab = "Time (hrs)")

> sessionInfo()
R version 3.3.3 (2017-03-06)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.2 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] bindrcpp_0.1         vpc_0.1.1            magrittr_1.5         classInt_0.1-23      dplyr_0.5.0.9000     purrr_0.2.2         
 [7] readr_1.1.0          tidyr_0.6.1.9000     tibble_1.2           ggplot2_2.2.1        tidyverse_1.1.1.9000

loaded via a namespace (and not attached):
 [1] reshape2_1.4.2      splines_3.3.3       haven_1.0.0.9000    lattice_0.20-34     colorspace_1.3-2    htmltools_0.3.5     yaml_2.1.14        
 [8] survival_2.41-2     rlang_0.0.0.9004    e1071_1.6-8         foreign_0.8-67      glue_0.0.0.9000     DBI_0.6             modelr_0.1.0       
[15] readxl_0.1.1        plyr_1.8.4          bindr_0.1           stringr_1.2.0       munsell_0.4.3       gtable_0.2.0        rvest_0.3.2        
[22] psych_1.6.12        evaluate_0.10       labeling_0.3        knitr_1.15.1        forcats_0.2.0       parallel_3.3.3      class_7.3-14       
[29] broom_0.4.2         Rcpp_0.12.10        scales_0.4.1        backports_1.0.5     jsonlite_1.3        mnormt_1.5-5        hms_0.3            
[36] digest_0.6.12       stringi_1.1.3       grid_3.3.3          rprojroot_1.2       tools_3.3.3         lazyeval_0.2.0.9000 MASS_7.3-45        
[43] Matrix_1.2-8        xml2_1.1.1          lubridate_1.6.0     assertthat_0.1      rmarkdown_1.3       httr_1.2.1          R6_2.2.0           
[50] nlme_3.1-131 

[ronkeizer]

if you use "classic" data.frames for sim/obs instead of tbl_df it works fine. Will put a conversion in place to do this automatically.

[Eliford]

Thank you Ron but the solution still does not work for me. I converted from tbl_df to data.frame as below. But still got the error! Error in grep("row", facet) : object 'facet' not found

obs<-obs%>%mutate(doseg=ifelse(dose<4.5,1,2))%>%as.data.frame()
sim<-sim%>%mutate(doseg=ifelse(dose<4.5,1,2))%>%as.data.frame()

[ronkeizer]

i already implemented the auto-conversion, so if you just re-install the vpc library from GitHub the code you sent earlier should work fine even without converting to data.frame yourself.

[Eliford]

Dear Ron, Thank you for working on this issue. Unfortunately I still get the same error message as before. Below is the session info when I run the same code as I posted before.

> sessionInfo()
R version 3.3.3 (2017-03-06)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.2 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] bindrcpp_0.1     vpc_0.1.1        magrittr_1.5     classInt_0.1-23  dplyr_0.5.0.9001

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.10        bindr_0.1           knitr_1.15.1        splines_3.3.3       MASS_7.3-45         munsell_0.4.3       colorspace_1.3-2   
 [8] lattice_0.20-34     R6_2.2.0            rlang_0.0.0.9006    plyr_1.8.4          tools_3.3.3         grid_3.3.3          gtable_0.2.0       
[15] e1071_1.6-8         class_7.3-14        yaml_2.1.14         survival_2.41-2     lazyeval_0.2.0.9000 assertthat_0.1      tibble_1.3.0       
[22] Matrix_1.2-8        ggplot2_2.2.1       glue_0.0.0.9000     scales_0.4.1    

[ronkeizer]

ah i see it now, issue with "facet" variable in the plot_vpc(). I had a global "facet" variable set in my environment so I didn't see this bug. Fixed now, please re-install.

[Eliford]

Thank you Ron, i think you haven't pushed the new commit.

[ronkeizer]

done now :)