run1 <- PM_load(1, path = "Runs")8 Examine results
Recall that in the previous exercises, we loaded our run results with the following.
run1 is now a PM_result R6 object that contains all the results from the NPAG run. Let’s dig into it in more detail.
8.1 Objects within a PM_result
8.1.1 Data
8.1.1.1 Data plots
You can plot the original data using R6 with various options. Type ?plot.PM_data in the R console for help.
run1$data$plot() Here’s some nice grouping and coloring. If you use grouping, provide names in group_names and colors/symbols for each group as arguments to marker. Colors of lines will be set to the same as for marker to avoid a mess of colors.
run1$data$plot(group = "gender", group_names = c("Female", "Male"),
marker = list(color = c("red","blue"), symbol = c("circle","triangle-up"))) The next one splits the data by subject but will not render easily in this document. You can interact with it in your own script or Quarto document.
run1$data$plot(overlay = FALSE, xlim = c(119, 145))The following are the older S3 method with plot(...) for the first two examples. You can use R6 or S3 for any Pmetrics object. We will focus on R6 as the more modern way.
plot(run1$data)8.1.1.2 Data summaries
Here’s a summary of the original data file; ?summary.PM_data for help.
run1$data$summary()8.1.2 Observed vs. Predicted (OP)
8.1.2.1 OP plots
You can plot observed vs. predicted data. Type ?plot.PMop in the R console for help or check out plot.PM_op.
run1$op$plot() Below, you can see how to change the prediction type, regression line and color, and plotting symbol and color. Most plots in Pmetrics have marker and line list arguments to customize the appearance of the plotting symbols and lines.
run1$op$plot(pred.type = "pop")
run1$op$plot(line = list(lm = FALSE, loess = list(color = "red")), marker = list(symbol = 3, color = "green"))8.1.2.2 OP summaries
Get a summary with bias and imprecision of the population predictions; ?summary.PM_op in the console or summary.PM_op for more information. The default pred.type is “post”, which are the predictions based on the full Bayesian posterior distribution of parameter values for an individual. The predictions are the weighted mean or median (controlled by the icen argument in all Pmetrics functions) of the predictions from each support point in the individual’s Bayesian posterior, joint probability parameter distribution, i.e. the support points in the posterior. Here we will look at the population predictions based on the means of the parameter values.
run1$op$summary(pred.type = "pop", icen = "mean")The nicely formatted table below will open in your Viewer panel giving you the summary statistics for non-missing observations and corresponding predictions.
Notes at the bottom of the table tell you the basis for the predictions, which can be changed with arguments to the $summary() method, which calls summary.PM_op. For example, aince neither “pop” nor “mean” are the defaults, we have to specify them. We can also use the S3 method for the same summary below.
summary(run1$op, pred.type = "pop", icen = "mean") # S3 methodAs with any object in Pmetrics, the PM_op R6 object has fields and methods. As a reminder of what we’ve covered in several previous chapters, R6 fields contain data and methods are functions that act on the data. All fields and methods are accessed with $ added to the name of the object. For example, in a PM_op object and most other Pmetrics objects, the raw data can be accessed via the $data field.
run1$op$dataHaving the raw data provides a data frame compatible with any R function, particularly those from Tidyverse, including piping (%>%). If you’re not familiar with the pipe operator, see this article. We strongly recommend Tidvyerse for data manipulation and visualization in R, and Pmetrics works well with it. See R for Data Science for more information.
This allows pre-processing in ways more flexible than the default plot method.
library(tidyverse)
run1$op$data %>% plot()
run1$op$data %>%
filter(pred > 5) %>%
filter(pred < 10) %>%
summary()See a header with the first 10 rows of the op object:
head(run1$op$data, 10)8.1.3 Final population joint density
8.1.3.1 Final plots
Plot the final population joint density. Type ?plot.PM_final in the R console for help.
run1$final$plot()
# alternate methods with same result
run1$final$data %>% plot() # tidyverse style
plot(run1$final) # S3 method
plot(run1$final$data) # S3 methodYou can add a kernel density curve.
run1$final$plot(line = list(color = "red"))It is possible to make a bivariate plot. Plotting formulae in R are of the form y~x
run1$final$plot(ke ~ v,
marker = list(color = "red", symbol = "diamond"))8.1.3.2 Final summaries
See a summary with confidence intervals around the medians and the Median Absolute Weighted Difference (MAWD); ?summary.PM_final for help and further explanation.
run1$final$summary()The original final object can be accessed with the $data field.
run1$final$data
names(run1$final$data)See the population points:
run1$final$popPoints
#> # A tibble: 19 × 5
#> ka ke v tlag1 prob
#> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0.899 0.0340 74.4 2.01 0.101
#> 2 0.898 0.0433 108. 2.09 0.0500
#> 3 0.897 0.0608 32.5 3.14 0.0500
#> 4 0.899 0.0564 120. 0.678 0.0502
#> 5 0.363 0.0749 68.4 0.869 0.0500
#> 6 0.805 0.0755 60.8 1.09 0.0500
#> 7 0.899 0.0278 63.7 1.75 0.0506
#> 8 0.690 0.0371 85.1 0.644 0.0504
#> 9 0.102 0.0642 119. 0.00843 0.0498
#> 10 0.717 0.0590 64.8 0.821 0.0497
#> 11 0.364 0.0278 58.9 0.458 0.0494
#> 12 0.549 0.0450 114. 0.316 0.0521
#> 13 0.694 0.0456 97.9 0.853 0.0497
#> 14 0.748 0.0291 95.9 1.39 0.0502
#> 15 0.866 0.0331 74.4 1.11 0.0484
#> 16 0.747 0.0397 120. 0.00315 0.0485
#> 17 0.437 0.0698 43.2 3.99 0.0568
#> 18 0.727 0.0578 53.3 3.04 0.0149
#> 19 0.725 0.0578 53.3 3.04 0.0783Or drill down one more level into the raw data. The $final$popPoints field above pulls its information from the $final$data$popPoints field below, but is provided as a simpler way to access the population points.
run1$final$data$popPointsSee the population mean parameter values:
run1$final$popMean
#> # A tibble: 1 × 4
#> ka ke v tlag1
#> <dbl> <dbl> <dbl> <dbl>
#> 1 0.694 0.0487 79.0 1.488.1.4 Cycle information
8.1.4.1 Cycle plots
Plot cycle information; type ?plot.PM_cycle in the R console for help.
run1$cycle$plot()8.1.4.2 Cycle summaries
Summarize the cycle information; ?summary.PM_cycle for help.
run1$cycle$summary()
# alternatives
run1$cycle$data %>% summary() # tidyverse style
summary(run1$cycle) # S3 method
summary(run1$cycle$data) # S3 method#>
#> ── Cycle Summary ──
#>
#> Cycle number: 100 of 100
#> Status: Stop: MaxCycles
#> Log-likelihood: 440.985
#> AIC:: 448.985
#> BIC: 452.968
#> Outeq 1: Gamma = 2.39
#>
#> ── Normalized parameter values:
#> # A tibble: 3 × 6
#> cycle ka ke v tlag1 stat
#> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 100 1.14 0.978 0.975 0.931 mean
#> 2 100 1.10 0.622 1.83 1.41 sd
#> 3 100 1.13 1 0.92 0.865 median8.1.5 Covariate information
8.1.5.1 Covariate plots
Plot covariate information; type ?plot.PM_cov in the R console for help. Recall that plotting formulae in R are of the form y~x.
run1$cov$plot(v ~ wt)You can also use the raw data in the $data field to make plots with other R functions, including those from Tidyverse. Here we will filter the data for age > 25 and then plot volume vs. weight.
run1$cov$data %>%
filter(age > 25) %>%
plot(v ~ wt)Change the formatting of the plot with the line and marker list arguments.
run1$cov$plot(ke ~ age, line = list(loess = FALSE, lm = list(color = "red")),
marker = list(symbol = 3))Another plot with mean Bayesian posterior parameter and covariate values…remember the icen argument?
run1$cov$plot(v ~ wt, icen = "mean")When time is the x variable, the y variable is aggregated by subject. In R plot formulae, calculations on the fly can be included using the I() function.
run1$cov$plot(I(v * wt) ~ time)The above plots the product of volume and weight vs. time for each subject.
8.1.5.2 Covariate summaries
Just as you have the option to plot by mean covariate and posterior parameter values, you can summarize with means; ?summary.PM_cov for help.
run1$cov$summary(icen = "mean")
#> # A tibble: 20 × 11
#> id ka ke v tlag1 africa age gender height wt icen
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 1 0.370 0.0278 59.0 0.473 1 21 1 160 46.7 mean
#> 2 2 0.739 0.0399 120. 0.0155 1 30 1 174 66.5 mean
#> 3 3 0.898 0.0433 108. 2.09 1 24 0 164 46.7 mean
#> 4 4 0.899 0.0564 120. 0.678 1 25 1 165 50.8 mean
#> 5 5 0.104 0.0642 119. 0.0105 1 22 1 181 65.8 mean
#> 6 6 0.898 0.0340 74.4 2.01 1 23 1 177 65 mean
#> 7 7 0.706 0.0586 52.7 3.11 1 27 0 161 51.7 mean
#> 8 8 0.898 0.0340 74.4 2.01 1 22 1 163 51.2 mean
#> 9 9 0.692 0.0456 98.1 0.845 1 23 1 174 55 mean
#> 10 10 0.717 0.0590 64.9 0.821 1 32 1 163 52.1 mean
#> 11 11 0.363 0.0749 68.4 0.869 1 34 1 165 56.5 mean
#> 12 12 0.749 0.0292 95.5 1.39 1 54 0 160 47.9 mean
#> 13 13 0.452 0.0692 43.7 3.94 1 24 1 180 60.5 mean
#> 14 14 0.552 0.0449 114. 0.315 1 26 1 174 59.2 mean
#> 15 15 0.864 0.0330 74.8 1.13 1 19 0 150 43 mean
#> 16 16 0.690 0.0371 85.2 0.646 1 25 1 173 64.4 mean
#> 17 17 0.805 0.0755 60.8 1.09 1 23 1 170 54.8 mean
#> 18 18 0.899 0.0278 63.7 1.74 1 20 0 164 44.3 mean
#> 19 19 0.897 0.0608 32.5 3.14 1 36 1 168 50 mean
#> 20 20 0.691 0.0592 52.1 3.16 1 31 1 170 59 meanThe run1$cov object looks like a data frame, but is really an R6 PM_cov object whose $print() method returns a data frame. You can access the raw data with the $data field, which is truly a data frame.
run1$cov$data[, 1:3] # for example
run1$cov$data %>% filter(gender == 1) %>% summary() # using TidyverseThis returns a data frame with the means of the individual’s covariates over the observation period and the mean of the Bayesian posterior parameter values.
When trying to find covariate-parameter relationships, it is convenient to examine at all possible covariate-parameter relationships by multiple linear regression with forward and backward elimination - type ?PM_step in the R console for help or see PM_step.
run1$step()
#> age gender height wt
#> ka NA NA NA NA
#> ke NA 0.047 0.076 0.105
#> v NA NA 0.181 0.116
#> tlag1 NA 0.162 NA NA
#> africa NA NA NA NA…or on the cov object directly.
run1$cov$step()icen works here too.
run1$step(icen = "mean")If you wish to see P values for forward elimination only:
run1$step(direction = "forward")
#> age gender height wt
#> ka 0.841 0.519 0.257 0.267
#> ke 0.920 0.127 0.236 0.281
#> v 0.784 0.694 0.084 0.149
#> tlag1 0.684 0.287 0.238 0.624
#> africa NA NA NA NA