Analysis with Julia

After you have Julia installed (see Julia Quick Start), you might want to analyze some data using Julia. Here’s a collection of examples that can help you do that. They are sorted based on a common workflow:

1. Import data
2. Plot data Many 2D and 3D plot examples have been collected in Beautiful Makie, and it’s useful to browse through this collection to understand ways to plot your data.
3. Fit data
4. Determine significance

Importing data

Here are some strategies to import and save data.

Importing data from a csv file

Use the CSV.jl and DataFrames.jl packages:

using CSV, DataFrames

n = 10 #Let's generate 10 random numbers
x = sort((rand( -10:100, n)))
y = 5/9 .* (x  .- 32) #A conversion from Fahrenheit to Celsius

data = DataFrame(F=x,C=y) #Turn it into a neat little data frame

CSV.write("data.csv", data) #write the data to a CSV
data_new = CSV.read("data.csv", DataFrame) #read the data from a CSV

x_values = data_new.F # Access columns by name
y_values = data_new.C

Importing a tiff image

Use TiffImages.jl for TIFF loading:

using TiffImages

img = TiffImages.load("image.tif")
height, width = size(img)
#Julia uses complicated types for images. I often convert to a simpler format.
img_float = Int16.(img) 

Importing data from a multidimensional tiff image

For stacks, time series, and multi-channel images:

using TiffImages

img_stack = TiffImages.load("stack.tif")

n_frames = size(img_stack, 3) # For 3D data (height × width × frames)

CHANNELS = 4 #As an example
channel_1 = img_stack[:, :, 1:CHANNELS:end]
channel_2 = img_stack[:, :, 2:CHANNELS:end]
channel_3 = img_stack[:, :, 3:CHANNELS:end]
channel_4 = img_stack[:, :, 4:CHANNELS:end]

Saving your Julia workspace with JLD2.jl

Save and load Julia variables and workspaces as .jdl2 files, like you can in Matlab using .mat files.

using JLD2

@save "workspace.jld2" # Save everything

# Save specific variables
@save "workspace.jld2" time_points measurements all_data

# Load variables back
@load "workspace.jld2" time_points measurements all_data

# Or load into a dictionary
data = load("results.jld2")
loaded_time = data["experiment_time"]

Importing data from an xml file

Use EzXML.jl for XML parsing. This is more difficult and probably requires browsing the XML file for a bit to understand it’s structure.

using EzXML

doc = parsexml(read("data.xml", String))

root = doc.root # Navigate to specific elements
for element in eachelement(root)
    value = element["attribute_name"]     # Extract attributes
    text = nodecontent(element)     # Extract text content
end

Fitting data

Use LsqFit.jl for curve fitting:

Exponential decay fit

using GLMakie, LsqFit

t_data = 0:9
y_data = [95.2, 74.1, 58.3, 46.2, 37.1, 30.2, 25.1, 21.5, 18.9, 17.1]

#= Exponential decay: y = A * exp(-t/τ) + baseline
We create a function we call exp_decay of that form
But this could be any function you want to fit to!=#
exp_decay(t, p) = p[1] .* exp.(-t ./ p[2]) .+ p[3]

p_guess = [100.0, 3.0, 10.0] # Initial parameter guesses

fit_result = curve_fit(exp_decay, t_data, y_data, p_guess)
fit_params = coef(fit_result)

#Fit with 100 subsamples
x_fit_range = range(t_data[1],t_data[end], 100)
y_fit = exp_decay(x_fit_range, fit_params)

#Plot
fig = Figure(size = (600, 400))
ax = Axis(fig[1, 1])
scatter!(ax, t_data, y_data, markersize = 10, color=:black)
lines!(ax, x_fit_range, y_fit, linewidth = 3, color=:red)

fig

Precompiling packages...
    520.8 ms  ✓ StaticArraysCore
    577.4 ms  ✓ Adapt
    881.6 ms  ✓ ProgressMeter
    888.8 ms  ✓ AxisArrays
    948.1 ms  ✓ FillArrays
   1302.3 ms  ✓ PDMats
    719.4 ms  ✓ Preferences
    832.5 ms  ✓ SimpleTraits
   2104.6 ms  ✓ IrrationalConstants
    645.2 ms  ✓ Adapt → AdaptSparseArraysExt
    543.4 ms  ✓ OffsetArrays → OffsetArraysAdaptExt
   1537.2 ms  ✓ Compat
   1889.3 ms  ✓ StructArrays
    574.5 ms  ✓ FillArrays → FillArraysStatisticsExt
    838.1 ms  ✓ FillArrays → FillArraysSparseArraysExt
    645.7 ms  ✓ PrecompileTools
    813.9 ms  ✓ FillArrays → FillArraysPDMatsExt
    894.9 ms  ✓ JLLWrappers
   4594.5 ms  ✓ BaseDirs
   1226.6 ms  ✓ StructArrays → StructArraysAdaptExt
   1298.9 ms  ✓ LogExpFunctions
   1955.2 ms  ✓ ComputePipeline
   1289.6 ms  ✓ Compat → CompatLinearAlgebraExt
   3357.9 ms  ✓ FileIO
   1686.4 ms  ✓ StructArrays → StructArraysSparseArraysExt
    532.2 ms  ✓ StructArrays → StructArraysLinearAlgebraExt
    689.5 ms  ✓ Graphite2_jll
    640.2 ms  ✓ Libmount_jll
    636.7 ms  ✓ EpollShim_jll
    694.3 ms  ✓ LLVMOpenMP_jll
    692.5 ms  ✓ Bzip2_jll
    704.4 ms  ✓ Rmath_jll
    649.2 ms  ✓ Xorg_libXau_jll
    745.0 ms  ✓ libpng_jll
    766.4 ms  ✓ libfdk_aac_jll
   3492.5 ms  ✓ ColorSchemes
    701.6 ms  ✓ Imath_jll
    673.9 ms  ✓ Giflib_jll
    699.2 ms  ✓ LAME_jll
   1175.1 ms  ✓ IntelOpenMP_jll
    690.4 ms  ✓ LERC_jll
    692.6 ms  ✓ EarCut_jll
    680.2 ms  ✓ CRlibm_jll
    762.4 ms  ✓ JpegTurbo_jll
    743.3 ms  ✓ XZ_jll
    681.3 ms  ✓ Ogg_jll
    754.2 ms  ✓ oneTBB_jll
    610.8 ms  ✓ Xorg_libXdmcp_jll
    846.8 ms  ✓ x265_jll
   7453.1 ms  ✓ StaticArrays
   1160.0 ms  ✓ x264_jll
   1083.5 ms  ✓ libaom_jll
    914.2 ms  ✓ Zstd_jll
   1083.9 ms  ✓ Expat_jll
    972.0 ms  ✓ Opus_jll
    980.7 ms  ✓ LZO_jll
    936.8 ms  ✓ Xorg_xtrans_jll
   1034.6 ms  ✓ Libiconv_jll
    998.3 ms  ✓ Libffi_jll
    980.2 ms  ✓ isoband_jll
   9963.1 ms  ✓ SIMD
   1224.4 ms  ✓ FFTW_jll
  10115.2 ms  ✓ Parsers
    663.9 ms  ✓ Libuuid_jll
    957.5 ms  ✓ OpenSpecFun_jll
    780.6 ms  ✓ FriBidi_jll
    810.4 ms  ✓ OpenBLASConsistentFPCSR_jll
    568.9 ms  ✓ LogExpFunctions → LogExpFunctionsInverseFunctionsExt
    723.0 ms  ✓ Pixman_jll
    790.6 ms  ✓ FreeType2_jll
    909.4 ms  ✓ Rmath
   1580.8 ms  ✓ QOI
    866.5 ms  ✓ CRlibm
   1109.3 ms  ✓ OpenEXR_jll
   2277.1 ms  ✓ DataStructures
   1015.6 ms  ✓ libsixel_jll
   2604.6 ms  ✓ ChainRulesCore
  17591.4 ms  ✓ Unitful
    765.4 ms  ✓ Xorg_libxcb_jll
   1003.3 ms  ✓ libvorbis_jll
    759.6 ms  ✓ StaticArrays → StaticArraysStatisticsExt
    750.9 ms  ✓ ConstructionBase → ConstructionBaseStaticArraysExt
    685.7 ms  ✓ Adapt → AdaptStaticArraysExt
   1603.5 ms  ✓ MKL_jll
    616.3 ms  ✓ Dbus_jll
   1016.0 ms  ✓ Libtiff_jll
    725.4 ms  ✓ Wayland_jll
    977.4 ms  ✓ GettextRuntime_jll
    680.2 ms  ✓ Isoband
   1792.1 ms  ✓ StructArrays → StructArraysStaticArraysExt
  16013.4 ms  ✓ ImageCore
   1825.7 ms  ✓ FreeType
   1697.8 ms  ✓ Fontconfig_jll
   1991.2 ms  ✓ JSON
    768.5 ms  ✓ SortingAlgorithms
   3086.7 ms  ✓ SpecialFunctions
   1902.8 ms  ✓ OpenEXR
   2101.1 ms  ✓ QuadGK
    554.1 ms  ✓ AbstractFFTs → AbstractFFTsChainRulesCoreExt
   1643.9 ms  ✓ ChainRulesCore → ChainRulesCoreSparseArraysExt
   3980.1 ms  ✓ IntervalArithmetic
    874.7 ms  ✓ Unitful → ConstructionBaseUnitfulExt
   1647.4 ms  ✓ LogExpFunctions → LogExpFunctionsChainRulesCoreExt
    845.9 ms  ✓ Unitful → InverseFunctionsUnitfulExt
   1839.5 ms  ✓ StaticArrays → StaticArraysChainRulesCoreExt
    906.3 ms  ✓ Xorg_libX11_jll
    927.6 ms  ✓ Unitful → PrintfExt
  10379.0 ms  ✓ PlotUtils
   1192.4 ms  ✓ Glib_jll
   1822.6 ms  ✓ ImageBase
    930.0 ms  ✓ ColorBrewer
   2818.7 ms  ✓ JpegTurbo
   2405.4 ms  ✓ Sixel
   4437.9 ms  ✓ PNGFiles
   1319.6 ms  ✓ HypergeometricFunctions
  11699.5 ms  ✓ Automa
   2165.2 ms  ✓ SpecialFunctions → SpecialFunctionsChainRulesCoreExt
  13229.0 ms  ✓ GeometryBasics
   1061.5 ms  ✓ IntervalArithmetic → IntervalArithmeticSparseArraysExt
   1085.1 ms  ✓ ColorVectorSpace → SpecialFunctionsExt
    748.8 ms  ✓ IntervalArithmetic → IntervalArithmeticIntervalSetsExt
   1012.4 ms  ✓ IntervalArithmetic → IntervalArithmeticLinearAlgebraExt
    758.4 ms  ✓ Xorg_libXext_jll
    780.7 ms  ✓ Xorg_libXfixes_jll
    793.0 ms  ✓ Xorg_libXrender_jll
   3970.9 ms  ✓ StatsBase
   7472.5 ms  ✓ FFTW
    761.7 ms  ✓ Xorg_libxkbfile_jll
    900.0 ms  ✓ Packing
   3215.4 ms  ✓ Interpolations
   1908.4 ms  ✓ ImageAxes
   1885.2 ms  ✓ ShaderAbstractions
    812.1 ms  ✓ Xorg_libXinerama_jll
    741.5 ms  ✓ Libglvnd_jll
   3098.4 ms  ✓ StatsFuns
   2872.5 ms  ✓ MeshIO
    751.7 ms  ✓ Xorg_libXi_jll
    759.3 ms  ✓ Xorg_libXrandr_jll
   2604.7 ms  ✓ FreeTypeAbstraction
    740.2 ms  ✓ Xorg_libXcursor_jll
    714.1 ms  ✓ Xorg_xkbcomp_jll
   1624.3 ms  ✓ Cairo_jll
    723.5 ms  ✓ StatsFuns → StatsFunsInverseFunctionsExt
   1171.8 ms  ✓ Interpolations → InterpolationsUnitfulExt
   1347.1 ms  ✓ ImageMetadata
   1358.8 ms  ✓ libwebp_jll
   4003.4 ms  ✓ ExactPredicates
   2216.3 ms  ✓ StatsFuns → StatsFunsChainRulesCoreExt
   1800.0 ms  ✓ Xorg_xkeyboard_config_jll
   6934.0 ms  ✓ GridLayoutBase
   1618.9 ms  ✓ HarfBuzz_jll
    636.5 ms  ✓ xkbcommon_jll
   1097.6 ms  ✓ libass_jll
   2562.4 ms  ✓ Netpbm
    906.2 ms  ✓ Pango_jll
   2110.7 ms  ✓ WebP
    581.2 ms  ✓ libdecor_jll
    634.0 ms  ✓ GLFW_jll
   1866.7 ms  ✓ FFMPEG_jll
   5452.8 ms  ✓ MathTeXEngine
   5658.4 ms  ✓ Distributions
   1119.6 ms  ✓ GLFW
   4198.8 ms  ✓ DelaunayTriangulation
   1138.3 ms  ✓ Distributions → DistributionsTestExt
   1792.5 ms  ✓ Distributions → DistributionsChainRulesCoreExt
   1469.0 ms  ✓ KernelDensity
  30170.9 ms  ✓ TiffImages
    911.4 ms  ✓ ImageIO
 116142.0 ms  ✓ Makie
  50795.5 ms  ✓ GLMakie
  170 dependencies successfully precompiled in 226 seconds. 115 already precompiled.
Precompiling packages...
   1684.4 ms  ✓ QuartoNotebookWorkerLaTeXStringsExt (serial)
  1 dependency successfully precompiled in 2 seconds
Precompiling packages...
   1136.4 ms  ✓ QuartoNotebookWorkerJSONExt (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
   1047.7 ms  ✓ QuartoNotebookWorkerTablesExt (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
   5415.2 ms  ✓ QuartoNotebookWorkerMakieExt (serial)
  1 dependency successfully precompiled in 5 seconds
Precompiling packages...
    478.9 ms  ✓ DiffResults
    563.0 ms  ✓ CommonSubexpressions
    650.3 ms  ✓ ADTypes
    662.9 ms  ✓ DiffRules
   1039.6 ms  ✓ Setfield
    466.6 ms  ✓ ADTypes → ADTypesConstructionBaseExt
   1328.9 ms  ✓ ArrayInterface
    444.5 ms  ✓ ArrayInterface → ArrayInterfaceStaticArraysCoreExt
    539.7 ms  ✓ ArrayInterface → ArrayInterfaceSparseArraysExt
   1536.7 ms  ✓ DifferentiationInterface
    577.7 ms  ✓ FiniteDiff
    572.8 ms  ✓ DifferentiationInterface → DifferentiationInterfaceSparseArraysExt
    514.5 ms  ✓ DifferentiationInterface → DifferentiationInterfaceFiniteDiffExt
    569.9 ms  ✓ FiniteDiff → FiniteDiffSparseArraysExt
   3521.1 ms  ✓ ForwardDiff
   1418.5 ms  ✓ DifferentiationInterface → DifferentiationInterfaceForwardDiffExt
   1346.5 ms  ✓ NLSolversBase
   2650.5 ms  ✓ LsqFit
  18 dependencies successfully precompiled in 10 seconds. 58 already precompiled.
Precompiling packages...
    486.2 ms  ✓ IntervalArithmetic → IntervalArithmeticDiffRulesExt
  1 dependency successfully precompiled in 1 seconds. 43 already precompiled.
Precompiling packages...
    625.6 ms  ✓ Unitful → ForwardDiffExt
  1 dependency successfully precompiled in 1 seconds. 22 already precompiled.
Precompiling packages...
   1271.9 ms  ✓ IntervalArithmetic → IntervalArithmeticForwardDiffExt
  1 dependency successfully precompiled in 2 seconds. 48 already precompiled.
Precompiling packages...
   1048.2 ms  ✓ ForwardDiff → ForwardDiffStaticArraysExt
  1 dependency successfully precompiled in 1 seconds. 22 already precompiled.
Precompiling packages...
    437.8 ms  ✓ ADTypes → ADTypesChainRulesCoreExt
  1 dependency successfully precompiled in 1 seconds. 9 already precompiled.
Precompiling packages...
    441.3 ms  ✓ DifferentiationInterface → DifferentiationInterfaceChainRulesCoreExt
  1 dependency successfully precompiled in 1 seconds. 11 already precompiled.
Precompiling packages...
    534.4 ms  ✓ DifferentiationInterface → DifferentiationInterfaceStaticArraysExt
  1 dependency successfully precompiled in 1 seconds. 10 already precompiled.
Precompiling packages...
    431.3 ms  ✓ ArrayInterface → ArrayInterfaceChainRulesCoreExt
  1 dependency successfully precompiled in 1 seconds. 11 already precompiled.
Precompiling packages...
    527.4 ms  ✓ FiniteDiff → FiniteDiffStaticArraysExt
  1 dependency successfully precompiled in 1 seconds. 21 already precompiled.

Determining Significance

Significance is a genuinely helpful tool to help you think about your data critically. The greater scientific community demands the simplicity of presented p-values often in the form of dazzling and deceitful stars ***.

Here I focus on both the bread and butter significance tests, and will add very specific cases where I find significance calculations a little more sticky.

Determining significance between two sets of numbers

Use HypothesisTests.jl for statistical tests. Here we show a t-test unpaired data

using GLMakie, HypothesisTests

#Two independent groups
group1 = [120, 125, 130, 128, 135]
group2 = [100, 105, 110, 108, 115] #AKA Group W

#Perform unpaired t-test (Welch's t-test)
p_val = pvalue(UnequalVarianceTTest(group1, group2))

#a little funtion to get those dazzling stars
stars(p) = p < 0.001 ? "***" : p < 0.01 ? "**" : p < 0.05 ? "*" : "ns"

#Plot
fig = Figure(size = (600, 600))
ax = Axis(fig[1, 1], xticks = ([1, 2], ["Group 1", "Group W"]), ylabel = "Measurement")

#This thing (x -> 1)., is a micro-function that takes everything and sets it to 1 applied dot-wise
boxplot!(ax, (x -> 1).(group1) , group1, color = :lightblue, width = 0.4)
boxplot!(ax, (x -> 2).(group1) , group2, color = :lightcoral, width = 0.4)
scatter!(ax, (x -> 1).(group2) , group1, color = :blue)
scatter!(ax, (x -> 2).(group2) , group2, color = :red)

#Significance Bracket
y_max = max(maximum(group1), maximum(group1)) * 1.1
bracket_points = [(1, y_max - 2), (1, y_max), (2, y_max), (2, y_max - 2)]
lines!(ax, bracket_points, color = :black)
text!(ax, 1.5, y_max, text=stars(p_val), align=(:center, :bottom),fontsize=20)

#Adjust y-axis limits to ensure everything is visible
ylims!(ax, nothing, y_max * 1.05)

fig

Precompiling packages...
   2891.1 ms  ✓ Roots
   1731.0 ms  ✓ HypothesisTests
  2 dependencies successfully precompiled in 5 seconds. 61 already precompiled.
Precompiling packages...
    470.2 ms  ✓ Accessors → StructArraysExt
  1 dependency successfully precompiled in 1 seconds. 20 already precompiled.
Precompiling packages...
   1078.1 ms  ✓ SimpleBufferStream
    800.0 ms  ✓ BitFlags
   1084.4 ms  ✓ DelimitedFiles
   1090.9 ms  ✓ CodecZlib
   1119.6 ms  ✓ Measures
   1143.9 ms  ✓ URIs
   1678.6 ms  ✓ Unzip
    642.1 ms  ✓ Xorg_libICE_jll
    703.6 ms  ✓ Accessors → UnitfulExt
    782.2 ms  ✓ LoggingExtras
    762.9 ms  ✓ ExceptionUnwrapping
   1033.1 ms  ✓ ConcurrentUtilities
    753.4 ms  ✓ fzf_jll
   2115.1 ms  ✓ Crayons
    691.6 ms  ✓ mtdev_jll
    665.4 ms  ✓ libevdev_jll
    686.2 ms  ✓ eudev_jll
    845.2 ms  ✓ MbedTLS_jll
    665.6 ms  ✓ Xorg_xcb_util_jll
    706.6 ms  ✓ Accessors → StaticArraysExt
    760.2 ms  ✓ Vulkan_Loader_jll
    728.7 ms  ✓ FFMPEG
   1510.2 ms  ✓ RecipesBase
    835.7 ms  ✓ Xorg_libSM_jll
    808.8 ms  ✓ libinput_jll
    893.5 ms  ✓ JLFzf
    881.1 ms  ✓ Xorg_xcb_util_image_jll
   1907.8 ms  ✓ OpenSSL
    814.4 ms  ✓ Xorg_xcb_util_keysyms_jll
    875.2 ms  ✓ Xorg_xcb_util_renderutil_jll
   1590.2 ms  ✓ MbedTLS
    796.3 ms  ✓ Xorg_xcb_util_wm_jll
   4709.1 ms  ✓ Latexify
    811.6 ms  ✓ Xorg_xcb_util_cursor_jll
   3502.8 ms  ✓ PlotThemes
   4165.7 ms  ✓ MarchingCubes
    901.3 ms  ✓ Latexify → SparseArraysExt
   1525.5 ms  ✓ Qt6Base_jll
   1693.7 ms  ✓ UnitfulLatexify
   3936.3 ms  ✓ RecipesPipeline
   2311.4 ms  ✓ Qt6ShaderTools_jll
   2471.5 ms  ✓ GR_jll
   1913.9 ms  ✓ Qt6Declarative_jll
    533.1 ms  ✓ Qt6Wayland_jll
  12490.5 ms  ✓ HTTP
   3532.4 ms  ✓ GR
  36014.6 ms  ✓ UnicodePlots
   1493.8 ms  ✓ UnicodePlots → UnitfulExt
  36381.7 ms  ✓ Plots
   2185.1 ms  ✓ Plots → UnitfulExt
  50 dependencies successfully precompiled in 61 seconds. 154 already precompiled.
Precompiling packages...
    552.0 ms  ✓ Roots → RootsForwardDiffExt
  1 dependency successfully precompiled in 1 seconds. 31 already precompiled.
Precompiling packages...
    484.2 ms  ✓ Roots → RootsChainRulesCoreExt
  1 dependency successfully precompiled in 1 seconds. 19 already precompiled.

Determining if a slope is significantly different from 0

When fitting a line, we test if the slope could be statistically the same as zero just by chance. The null hypothesis is that the slope is 0. So we calculate how many standard errors the slope is from zero? Pop-quiz then, in the inverse case: How do you prove the slope is significantly similar to 0?

using GLMakie, GLM

x_data = 1:10
y_data = [2, 2, 3, 3, 3, 3, 4, 4, 4, 5]
X = hcat((x->1).(x_data), x_data) 

fit = lm(X, y_data)
coeftable(fit) #the coeftable of fit has p-values in it's 4 column
p_val = coeftable(fit).cols[4][2]

x_fit = range(x_data[1]-1,x_data[end]+1,100)
X_fit = hcat((x->1).(x_fit), x_fit)

#Predict the fit line and the 95% conf. int.
preds = predict(fit, X_fit, interval = :confidence)

#Plot
fig = Figure(size=(600, 600))
ax = Axis(fig[1, 1], xlabel="x", ylabel="y")

scatter!(ax, x_data, y_data, color = :black)
lines!(ax, x_fit, preds.prediction, color = :blue, linewidth=2)
band!(ax, x_fit, preds.lower, preds.upper, color=(:skyblue, 0.2))

#p-value text
p_text = p_val < 0.0001 ? "p < 0.0001" : "p = $(round(p_val, digits=4))"
text!(ax, 0.05, 0.95, text=p_text, space=:relative, align=(:left, :top), fontsize=16)

fig

Precompiling packages...
    417.9 ms  ✓ ShiftedArrays
   1337.8 ms  ✓ StatsModels
   1466.9 ms  ✓ GLM
  3 dependencies successfully precompiled in 3 seconds. 53 already precompiled.

Determining if a slope is significantly 0

Ok so here is the tricky bit. What if it looks like there really is no relation, and we want to see if that is significant? Well… it’s hard. One thing you can do to make it easy is define a small range around zero where you would consider the slope to be functionally zero for your specific problem. Then you can test the hypothesis that it lies inside or outside of that range.

using GLM, Distributions

x_data = 1:20
y_data = [10.1, 9.9, 10.0, 10.2, 9.8, 10.1, 9.9, 10.0, 9.9, 10.1,
          9.8, 10.2, 10.0, 9.9, 10.1, 10.0, 9.8, 10.2, 9.9, 10.0]

#Let any slope between -0.1 and 0.1 be definitionally zero.
delta = 0.1

X = hcat(ones(length(x_data)), x_data)
fit = lm(X, y_data)

slope = coeftable(fit).cols[1][2]   
std_err = coeftable(fit).cols[2][2]
dof = dof_residual(fit)   # Degrees of freedom

#Perform the Two One-Sided T-Tests (requires Distributions.jl) ---
t_dist = TDist(dof) # Create a t-distribution object for our fit

#Is the slope significantly GREATER than the lower bound?
t_lower = (slope - (-delta)) / std_err
p_lower = 1 - cdf(t_dist, t_lower)

# Is the slope significantly LESS than the upper bound?
t_upper = (slope - delta) / std_err
p_upper = cdf(t_dist, t_upper)

#The final TOST p-value is the larger of the two
tost_p_value = max(p_lower, p_upper)

println("Equivalence Test Results (bounds = ±$delta):")
println("Slope: $(round(slope, digits=4)), SE: $(round(std_err, digits=4))")
println("TOST p-value: $(round(tost_p_value, digits=4))")

if tost_p_value < 0.05
    println("Conclusion: The slope is statistically equivalent to zero within delta=$delta.")
else
    println("Conclusion: We cannot conclude the slope is equivalent to zero within delta=$delta.")
end

Determining significance between two slopes

ANCOVA tests whether regression lines from different groups have different slope. This is a complex idea that’s worth spelling out here.

ANOVA tests if the means of two or more groups are statistically different.

The model is: \(Y_{ij} = \mu_i + \epsilon_{ij}\)

• \(Y_{ij}\) is the measured value for the j-th data point in the i-th group.
• \(\mu_i\) is the mean of group i.
• \(\epsilon_{ij}\) is the difference of a data point and the mean.

The null hypothesis (\(H_0\)) is that all group means are equal: \(H_0: \mu_1 = \mu_2 = \dots = \mu_k\)

This is tested statistically using the F-statistic:

\(F = \frac{\text{Variance between groups}}{\text{Variance within groups}}\)

A large F-statistic indicates the variation between groups (the signal) is much larger than the random variation within groups (the noise). A p-value then determines if this F-statistic is large enough to reject the null hypothesis.

ANCOVA extends this by testing if the group means are different after controlling for a continuous variable, the Covariate.

The model becomes: \(Y_{ij} = \mu_i + \beta(X_{ij} - \bar{X}) + \epsilon_{ij}\)

• \(\mu_i\) represents the adjusted mean of group i.
• \(\beta\) is the coefficient quantifying the covariate’s effect.
• \(X_{ij}\) is the covariate score for the j-th person in the i-th group.
• \(\bar{X}\) is the average of all covariate scores across all groups.

ANCOVA uses the same F-test logic to test two primary hypotheses:

• Null Hypothesis (\(H_0\)): All adjusted group means are equal. \(H_0: \mu_1 = \mu_2 = \dots = \mu_k\)

• Null Hypothesis (\(H_0\)): The regression coefficient for the covariate is zero (i.e., there is no linear relationship). \(H_0: \beta = 0\)

In code, this ANCOVA test is:

using GLM, DataFrames, GLMakie, StatsModels

#Define the data for each group separately for clarity
A_X = [1.538, 6.994, 1.019, 6.603, 7.603, 4.914, 8.773, 4.021, 7.186, 6.695]
A_Y = [9.424, 17.639, 6.949, 14.459, 21.308, 14.703, 25.158, 12.910, 23.893, 19.834]

B_X = [3.746, 5.175, 0.720, 9.524, 4.982, 6.769, 7.224, 9.479, 1.495, 1.249]
B_Y = [14.650, 19.031, 9.555, 29.480, 20.433, 23.441, 22.306, 25.264, 6.952, 12.501]

C_X = [0.560, 4.536, 8.414, 5.485, 2.953, 8.107, 2.994, 7.255, 5.984, 3.541]
C_Y = [3.650, 22.332, 41.558, 26.173, 15.339, 38.155, 16.305, 25.017, 27.204, 18.539]

#Combine the separate group data into a single DataFrame for analysis
data = vcat(
    DataFrame(X = A_X, Y = A_Y, Group = "A"),
    DataFrame(X = B_X, Y = B_Y, Group = "B"),
    DataFrame(X = C_X, Y = C_Y, Group = "C")
)

#Little function to get stars
stars(p) = p < 0.001 ? "***" : p < 0.01 ? "**" : p < 0.05 ? "*" : "ns"

# Function to reorder groups to help with multiple pairs
function reorder_groups(df, ref_group)
    ref_rows = filter(row -> row.Group == ref_group, df)
    other_rows = filter(row -> row.Group != ref_group, df)
    return vcat(ref_rows, other_rows)
end

# Fit models with different reference groups
data_A_ref = reorder_groups(data, "A")
ancova_A = lm(@formula(Y ~ 1 + X + Group + X & Group), data_A_ref)
ct_A = coeftable(ancova_A)

data_B_ref = reorder_groups(data, "B")
ancova_B = lm(@formula(Y ~ 1 + X + Group + X & Group), data_B_ref)
ct_B = coeftable(ancova_B)

get_p(model, term) = coeftable(model).cols[4][findfirst(==(term), coeftable(model).rownms)]
p_AB = get_p(ancova_A, "X & Group: B")
p_AC = get_p(ancova_A, "X & Group: C")
p_BC = get_p(ancova_B, "X & Group: C")

fig = Figure(size = (600, 600))
ax = Axis(fig[1, 1])

colors = Dict("A" => :blue, "B" => :orange, "C" => :green)
for group in ["A", "B", "C"]
    subset = filter(row -> row.Group == group, data)
    scatter!(ax, subset.X, subset.Y, color = colors[group])
    
    x_range = range(extrema(subset.X)..., length=100)
    predict_df = DataFrame(X = x_range, Group = group)
    y_pred = predict(ancova_A, predict_df)
    y_pred_float = [y for y in y_pred]

    lines!(ax, x_range, y_pred_float, color = colors[group], linewidth = 2.5)
    text!(ax, x_range[end], y_pred_float[end], text = " $group", align = (:left, :center), color = colors[group])
end

legend_text = """
Slope Significance
A vs B: $(stars(p_AB))
A vs C: $(stars(p_AC))
B vs C: $(stars(p_BC))
"""

# Create a container for the legend
legend_layout = GridLayout(tellwidth = false, tellheight = false, halign = 0.98, valign = 0.05)
fig[1, 1] = legend_layout

Box(legend_layout[1, 1], color = (:white, 0.85), strokecolor = :black, strokewidth = 1)
Label(legend_layout[1, 1], legend_text, padding = (10, 10, 8, 8))

fig

Precompiling packages...
    456.1 ms  ✓ InvertedIndices
    523.4 ms  ✓ PooledArrays
    695.4 ms  ✓ InlineStrings
   1313.1 ms  ✓ StringManipulation
   1497.3 ms  ✓ SentinelArrays
  11195.6 ms  ✓ PrettyTables
  29124.9 ms  ✓ DataFrames
  7 dependencies successfully precompiled in 42 seconds. 28 already precompiled.
Precompiling packages...
    468.3 ms  ✓ InlineStrings → ParsersExt
  1 dependency successfully precompiled in 1 seconds. 9 already precompiled.