NIBR data analysis

How-to
21

Regarding the log ticks:
Currently, Makie calculates log ticks the same way that it calculates normal ticks, just inside the log transformed range of values. That’s why you can get ticks like 10^0.5 as 0.5 is usually a nice tick value to have. If, however, you want only integer values, here’s a simple way to get that, which might at some point be a built-in option in Makie:

# make a new placeholder type
struct IntegerTicks end

# integer ticks are just all integers within the axis limits
# (this should work ok if there's only a couple orders of magnitude, otherwise one would want to reduce the number of ticks)
MakieLayout.get_tickvalues(::IntegerTicks, vmin, vmax) = ceil(Int, vmin) : floor(Int, vmax)

f = Figure(resolution = (900, 600))

for (i, x) in enumerate([0.01, 0.001, 0.0001])
    lines(
        f[1, i],
        range(x, 9.99, length = 200),
        axis = (
            yscale = log10,
            yticks = LogTicks(IntegerTicks()), # the LogTicks wrapper takes care of the log transform before finding the IntegerTicks
            ytickformat = xs -> string.(xs), # this should usually format the resulting values as normal decimal numbers without trailing zeros 
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            yminorticksvisible = true,
        )
    )
end

f

grafik
grafik900×600 34.4 KB

22

Thank you very much @vijay @storopoli I was able to generate a quality graph finally using your code. I am touched by the help I am receiving. I was also able to control the marker, line sizes, marker color. The only thing I could not successfully change was the linecolor (see code below). Markercolor works, but not linecolor, not sure why.

ot = observations_vs_time(pop_nca, 
                                axis = (xlabel = "Time (hours)", 
                                        ylabel = "Concentration (μg/mL)",
                                        yscale=log10, ytickformat=x -> string.(round.(x; digits=1)), 
                                        yticks = [0.1, 1, 10],
                                        ygridwidth = 3, 
                                        yminorticksvisible = true,
                                        yminorgridvisible = true,
                                        yminorticks = IntervalsBetween(10),
                                        xminorticksvisible = true,
                                        xminorgridvisible = true,
                                        xminorticks = IntervalsBetween(5),
                                        limits = (nothing, nothing, nothing, 20),
                                        spinewidth = 2),
                                        columns = 4, rows = 3,
                                        markersize = 8,
                                        markercolor = :green,
                                        linewidth = 2,
                                        linecolor = :black,
                                        facet = ( combinelabels = true,))
#
save("../code/practice/nibr/data/linear_pk_plots.pdf", ot)

image
image1598×1206 125 KB

23

@Michael Thank you. That was it. paginate = true was the culprit. The code that @vijay @storopoli provided did allow me to control the cosmetics; however it only produces graphs for the first 12 subjects (perhaps due to the columns=4, rows=3). There are a total of 50 subjects. If I increase columns/rows, the plots are squished. I am unsure as to first generating all 50 profiles; and then save() them. Any advice?
Bob

24

@julius Thank you. Unfortunately, I am unable to adapt the code to my data. My experience is more on the clinical side; and not on the software side. I barely manage. It’s my bad. The for loop etc is where I am lost, as to its application to my data. The code provided @vijay @storopoli seem to do the job, unless your code is more efficient.
Bob

25

The Vector{Figure} that paginate = true generates could be saved out to individual numbered PDFs using something like the following:

foreach(enumerate(pk_plots)) do (nth, plot)
    save("plot-$nth.pdf", plot)
end
26

Sorry for confusing you, the for loop was not related to your data but should just show how the tick finding works with three different y limits. The solution you were given before is fine if you know the y tickpositions beforehand. My solution is automatic, so it could also be useful. You only need this part

struct IntegerTicks end
MakieLayout.get_tickvalues(::IntegerTicks, vmin, vmax) = ceil(Int, vmin) : floor(Int, vmax)

And then pass the yticks and ytickformat keyword arguments as shown

        axis = (
            yscale = log10,
            yticks = LogTicks(IntegerTicks()), # the LogTicks wrapper takes care of the log transform before finding the IntegerTicks
            ytickformat = xs -> string.(xs), # this should usually format the resulting values as normal decimal numbers without trailing zeros
27

Aha! @Michael Thank you. I guess there is no equivalent of knitr so all graphs are in one file. Imagine 5 PDFs for every type of graph one produces. Any way to consolidate?
Bob

28

report(pk_plots; title = "my title") might be what you’re wanting in that case, which should generate a single PDF document with all the paginated plots.

29

Wow that would be great. How do I embed my graph with my specifications into this report?

30

Could you describe in more detail what you would like to embed in the generated PDF? Currently report has limited customization options, we will be expanding on it in future releases though.

31

@Michael sure. I would like to retain the plot attributes as below in the report. Can this done?

ot = observations_vs_time(pop_nca, 
                                axis = (xlabel = "Time (hours)", 
                                        ylabel = "Concentration (μg/mL)",
                                        yscale=log10, ytickformat=x -> string.(round.(x; digits=1)), 
                                        yticks = [0.1, 1, 10],
                                        ygridwidth = 3, 
                                        yminorticksvisible = true,
                                        yminorgridvisible = true,
                                        yminorticks = IntervalsBetween(10),
                                        xminorticksvisible = true,
                                        xminorgridvisible = true,
                                        xminorticks = IntervalsBetween(5),
                                        limits = (nothing, nothing, nothing, 20),
                                        spinewidth = 2),
                                        columns = 4, rows = 3,
                                        markersize = 8,
                                        markercolor = :green,
                                        linewidth = 2,
                                        linecolor = :black,
                                        facet = ( combinelabels = true,))
32

That should work fine. report uses the same output as save(..., ot) so anything that appears when you manually call save on a Figure should also appear in the PDF generated by report.

33

@Michael may be I am not communicating clearly here. I have 50 profiles, which have to spread out on multiple pages. Let us say 4 columns x 3 rows x 5 pages. Here is the code based on your advice I used:

#
pk_log_plots = observations_vs_time(pop_nca, 
                                axis = (xlabel = "Time (hours)", 
                                        ylabel = "Concentration (μg/mL)",
                                        yscale=log10, ytickformat=x -> string.(round.(x; digits=1)), 
                                        yticks = [0.1, 1, 10],
                                        ygridwidth = 3, 
                                        yminorticksvisible = true,
                                        yminorgridvisible = true,
                                        yminorticks = IntervalsBetween(10),
                                        xminorticksvisible = true,
                                        xminorgridvisible = true,
                                        xminorticks = IntervalsBetween(5),
                                        limits = (nothing, nothing, nothing, 20),
                                        spinewidth = 2),
                                        columns = 4, rows = 3,
                                        paginate = true,
                                        markersize = 8,
                                        markercolor = :green,
                                        linewidth = 2,
                                        linecolor = :black,
                                        facet = ( combinelabels = true,))
#
pk_log_plots[1]
pk_log_plots[2]
report(pk_log_plots; title = "pk_log_plots.PDF")

The pk_log_plots.PDF contains only the first page. Please note I used: paginate=true. Hence I am able to view pk_log_plots[1], [2]…[6].

34

hi @bobbrown
I need to check if this works in your version on juliahub.com but please try this

report(pk_log_plots, 
      title = "Individual fits", 
      output = "./foldertosaveplots",
      header = "Single Ascending Dose",
      footer = "Confidential")

check for the pdf in the working directory

35

As @vijay mentions above, output = "..." can be used to set where the report PDF appears, it defaults to a output folder in the current working directory. Is there anything in that folder (if it exists)?

36

@Michael @vijay We have a winner finally! I was able to output all the plots into this PDF file. Thank you all.

image
image3186×938 159 KB

Now I am getting greedy (and have nothing else to do). I also have another plot for average PK profiles: avg_pk_plots. I tried if I can concatenate the two sets of plots using this code. Of course, it did not work.

report([pk_log_plots; avg_pk_plots],
     #(pk_log_plots; avg_pk_plots)
     #(pk_log_plots, avg_pk_plots) 
      title = "Individual Semi-Log Exploratory Plots", 
      output = "./practice/nibr/data/plots",
      header = "Single Ascending Dose",
      footer = "Confidential")

I do want to do a pulse-check. Hope these discussions are useful for other users as well. You both seem to be pros. I am happy share my entire project work here, so others can learn as well. Not imposing anything. Just a thought.
Bob

37

@storopoli your code was immensely helpful. I still could not get linecolor = :red to work. Anything that you could see in my code that is erroneous?
B

38

It’ll be color = :red rather than linecolor = :red for observations_vs_time if I recall correctly.

39

That worked. I could not have a guessed it, given there is a markercolor argument. Thank you.

40

@bobbrown you need to do this

report(vcat(pk_log_plots, avg_pk_plots),
      title = "Summary and Individual Semi-Log Exploratory Plots", 
      output = "./practice/nibr/data/plots",
      header = "Single Ascending Dose",
      footer = "Confidential")

some other julia experts can weigh in here, but the pattern of concatenating two vectors in julia is that you need to split the vectors into individual elements using ... and then put them together. See the example below.

julia> a = [1,2,3]
3-element Vector{Int64}:
 1
 2
 3

julia> b = [4, 5, 6]
3-element Vector{Int64}:
 4
 5
 6

julia> [a, b]
2-element Vector{Vector{Int64}}:
 [1, 2, 3]
 [4, 5, 6]

julia> [a..., b...]
6-element Vector{Int64}:
 1
 2
 3
 4
 5
 6

#UPDATE: Based on Michael's response below
julia> vcat(a, b)
6-element Vector{Int64}:
 1
 2
 3
 4
 5
 6

The report function only takes a vector of plots and not a vector of vectors.

We are happy to have these discussions, Please keep them coming!

That would be awesome, please share your experiences. In the meanwhile, for those reading the thread, I summarized all the code required till now below.

using Pumas
using PumasUtilities
using Dates
using CairoMakie
using AlgebraOfGraphics
using DataFramesMeta
df = CSV.read("Single_Ascending_Dose_Dataset2.csv", DataFrame, missingstrings = ["NA"])
df = @rsubset df :TIME >= 0
@rtransform! df route = "ev"
##  map NCA data from CSV
ncapop = read_nca(df,
                id  = :ID, time = :NOMTIME, amt = :AMT, route = :route,
                observations = :LIDV, 
                group = [:DOSE])

#plot means               
avg_pk_plots = summary_observations_vs_time(ncapop, 
                                  color = :black, linewidth = 2, whiskerwidth = 8,
                                  paginate = true,
                                  limit = 1,
                                  axis = (xlabel = "Time (hours)", 
                                            ylabel = "Concentration (μg/mL)",
                                            yscale=log10, ytickformat=x -> string.(round.(x; digits=1)), 
                                            yticks = [0.1, 1, 10],
                                            ygridwidth = 3, 
                                            yminorgridcolor = :darkgrey,
                                            yminorticksvisible = true,
                                            yminorgridvisible = true,
                                            yminorticks = IntervalsBetween(10),
                                            xminorticksvisible = true,
                                            xminorgridvisible = true,
                                            xminorticks = IntervalsBetween(5),
                                            limits = (nothing, nothing, nothing, 30),
                                            spinewidth = 2),
                                    #facet = ( combinelabels = true,),
                                    figure = (resolution = (800,800),
                                                        fontsize = 22),)
                                  
pk_log_plots = observations_vs_time(ncapop, paginate = true, color = :black,
                                axis = (xlabel = "Time (hours)", 
                                        ylabel = "Concentration (μg/mL)",
                                        yscale=log10, ytickformat=x -> string.(round.(x; digits=1)), 
                                        yticks = [0.1, 1, 10],
                                        ygridwidth = 3, 
                                        yminorticksvisible = true,
                                        yminorgridvisible = true,
                                        yminorticks = IntervalsBetween(10),
                                        xminorticksvisible = true,
                                        xminorgridvisible = true,
                                        xminorticks = IntervalsBetween(5),
                                        limits = (nothing, nothing, nothing, 20),
                                        spinewidth = 2),
                                        columns = 4, rows = 4,
                                        facet = ( combinelabels = true,))

report([pk_log_plots..., avg_pk_plots...],
      title = "Summary and Individual Semi-Log Exploratory Plots", 
      output = "./practice/nibr/data/plots",
      header = "Single Ascending Dose",
      footer = "Confidential")

Hope this helps!