examine.qmd

---
title: "Examine results"
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  echo = TRUE,
  eval = FALSE
)
```

```{r}
#| label: setup
#| echo: false
#| message: false
#| eval: true
library(glue)
library(knitr)
library(Pmetrics)

r_help <- function(pkg, name, fn = NULL) {
    if(is.null(fn)) {fn <- name}
    glue::glue("[`{name}`](https://rdrr.io/pkg/{pkg}/man/{fn}.html)")
}

gh_help <- function(name) {

    glue::glue("[`{name}`](https://lapkb.github.io/Pmetrics/reference/{name}.html)")

}

pmetrics <- function(){
    knitr::asis_output("[Pmetrics]{style=\"color: #841010; font-family: 'Arial', Arial, sans-serif; font-weight: 900;\"}")
}


run1 <- NPex
```

Recall that in the previous exercises, we loaded our run results with the following.
```{r}
run1 <- PM_load(1, path = "Runs")
```

`run1` is now a `r gh_help("PM_result")` R6 object that contains all the results from the NPAG run. Let's dig into it in more detail.


## Objects within a PM_result

### Data

#### Data plots

You can plot the original data using R6 with various options. Type `?plot.PM_data` in the R console for help.

```{r}
#| eval: false
run1$data$plot() 
```

```{r}
#| eval: true
#| echo: false
plotly::plotly_build(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.

```{r}
#| eval: false
run1$data$plot(group = "gender", group_names = c("Female", "Male"), 
  marker = list(color = c("red","blue"), symbol = c("circle","triangle-up"))) 
```

```{r}
#| eval: true
#| echo: false
plotly::plotly_build(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.

```{r}
#| eval: false
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.

```{r}
#| eval: false
plot(run1$data)
```

If you would like to plot the observations and posterior predictions made at `idelta` intervals (see options for `$fit` method in `r gh_help("PM_model")` , you can do that by adding a `pred` line to a data plot. See `r gh_help("plot.PM_data")` for more information.

```{r}
#| eval: false
# suppress the segments joining observations and add posterior predictions with `pred = "post"`
run1$data$plot(line = list(join = FALSE, pred = "post"), xlim = c(119, 145)) 
```

Typically you would want to use `overlay = FALSE` with the above to avoid a mess of lines joining observations and predictions from different subjects, but it will not render easily in this document. You can interact with it in your own script or Quarto document.
```{r}
#| eval: true
#| echo: false
# suppress the segments joining observations and add posterior predictions with `pred = "post"`
plotly::plotly_build(run1$data$plot(line = list(join = FALSE, pred = "post"), xlim = c(119, 145))) 
```

### Observed vs. Predicted (OP)

#### OP plots

You can plot observed vs. predicted data. Type `?plot.PM_op` in the R console for help or check out `r gh_help("plot.PM_op")`.

```{r}
#| eval: false
run1$op$plot() 
```

```{r}
#| eval: true
#| echo: false
plotly::plotly_build(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.

```{r}
#| eval: false
run1$op$plot(pred.type = "pop")
run1$op$plot(line = list(lm = FALSE, loess = list(color = "red")), marker = list(symbol = 3, color = "green"))
```

```{r}
#| label: op_plots_options
#| eval: true
#| echo: false
plotly::plotly_build(run1$op$plot(pred.type = "pop"))
plotly::plotly_build(run1$op$plot(line = list(lm = FALSE, loess = list(color = "red")), marker = list(symbol = 3, color = "green")))
```

#### OP summaries

Get a summary with bias and imprecision of the population predictions; `?summary.PM_op` in the console or `r gh_help("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 `r 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.

```{r}
#| eval: false
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. 

![](Images/summaryOP.png)

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 `r gh_help("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.

```{r}
#| label: S3_op_summary
summary(run1$op, pred.type = "pop", icen = "mean") # S3 method
```

As with any object in `r 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 `r pmetrics()` objects, the **raw data** can be accessed via the `$data` field.

```{r}
run1$op$data
```

Having the raw data provides a data frame compatible with any R function, particularly those from [**Tidyverse**](https://www.tidyverse.org), including piping (`|>`). If you're not familiar with the pipe operator, see [this article](https://magrittr.tidyverse.org/). We strongly recommend Tidvyerse for data manipulation and visualization in R, and `r pmetrics()` works well with it. See [R for Data Science](https://r4ds.hadley.nz) for more information.

This allows pre-processing in ways more flexible than the default plot method.

```{r}
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:

```{r}
head(run1$op$data, 10)
```

### Final population joint density

#### Final plots

Plot the final population joint density. Type `?plot.PM_final` in the R console for help.

```{r}
run1$final$plot()

# alternate methods with same result
run1$final$data |> plot() # tidyverse style
plot(run1$final) # S3 method
plot(run1$final$data) # S3 method
```

```{r}
#| eval: true
#| echo: false
plotly::plotly_build(run1$final$plot())
```

You can add a kernel density curve.

```{r}
run1$final$plot(line = list(color = "red"))
```

```{r}
#| eval: true
#| echo: false
plotly::plotly_build(
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`

```{r}
run1$final$plot(ke ~ v,
  marker = list(color = "red", symbol = "diamond"))
```

```{r}
#| eval: true
#| echo: false
plotly::plotly_build(
run1$final$plot(ke ~ v,
  marker = list(color = "red", symbol = "diamond")))
```


#### 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.

```{r}
run1$final$summary()
```

The original final object can be accessed with the `$data` field.

```{r}
run1$final$data
names(run1$final$data)
```

See the population points:

```{r}
#| eval: true
run1$final$popPoints
```

Or 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.

```{r}
run1$final$data$popPoints
```

See the population mean parameter values:

```{r}
#| eval: true
run1$final$popMean
```


### Cycle information

#### Cycle plots

Plot cycle information; type `?plot.PM_cycle` in the R console for help.

```{r}
run1$cycle$plot()
```

```{r}
#| eval: true
#| echo: false
plotly::plotly_build(
run1$cycle$plot())
```

#### Cycle summaries


Summarize the cycle information; `?summary.PM_cycle` for help.

```{r}
run1$cycle$summary()

# alternatives
run1$cycle$data |> summary() # tidyverse style
summary(run1$cycle) # S3 method
summary(run1$cycle$data) # S3 method
```

```{r}
#| eval: true
#| echo: false
run1$cycle$summary()
```


### Covariate information

#### 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`.

```{r}
run1$cov$plot(v ~ wt)
```

```{r}
#| eval: true
#| echo: false
plotly::plotly_build(
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.

```{r}
run1$cov$data |>
  filter(age > 25) |>
  plot(v ~ wt)
```

```{r}
#| eval: false
#| echo: false

run1$cov$data |>
  dplyr::filter(age > 25) |>
  plot(v ~ wt)
```

Change the formatting of the plot with the `line` and `marker` list arguments.

```{r}
run1$cov$plot(ke ~ age, line = list(loess = FALSE, lm = list(color = "red")),
               marker = list(symbol = 3))
```

```{r}
#| eval: true
#| echo: false
plotly::plotly_build(
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?

```{r}
run1$cov$plot(v ~ wt, icen = "mean")
```

```{r}
#| eval: true
#| echo: false

plotly::plotly_build(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.

```{r}
run1$cov$plot(I(v * wt) ~ time)
```

```{r}
#| eval: true
#| echo: false
plotly::plotly_build(
run1$cov$plot(I(v * wt) ~ time))
```

The above plots the product of volume and weight vs. time for each subject.

#### 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.

```{r}
#| eval: true
run1$cov$summary(icen = "mean")
```


The `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.

```{r}
run1$cov$data[, 1:3] # for example
run1$cov$data |> filter(gender == 1) |> summary() # using Tidyverse
```


This 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 `r gh_help("PM_step")`.

```{r}
#| eval: true
run1$step()
```

...or on the cov object directly.

```{r}
run1$cov$step()
```

`icen` works here too.

```{r}
run1$step(icen = "mean")
```

If you wish to see P values for forward elimination only:

```{r}
#| eval: true
run1$step(direction = "forward")
```