Shortening lengthy code

New to Julia
1

Hello, I am trying to perform a linear regression for two separate layers using the same variables. Right now, I have written everything twice, once for layer one and once for layer two. Is there an easy way to simplify the code in order to have it execute the regression for both layers separately without having to write it twice?

using DataFrames, CSV, Statistics, GLMakie, DelimitedFiles, GLM

ResData = DataFrame(Layer = rand(1:2, 200), X = rand(200), Y = rand(200), Porosity = rand(200), Density = rand(200), Permeability = rand(200))

#renaming columns from sorted ResData to use later on 
LayerNum     = ResData[:, 1]
XCoord       = ResData[:, 2]
YCoord       = ResData[:, 3]
Porosity     = ResData[:, 4]
Density      = ResData[:, 5]
Permeability = ResData[:, 6]

#use mask to sort data into layers 1 and 2 
ResLayer1_mask = LayerNum .== 1 
ResLayer2_mask = LayerNum .== 2

l1x = XCoord[ResLayer1_mask]
l1y = YCoord[ResLayer1_mask]
l2x = XCoord[ResLayer2_mask]
l2y = YCoord[ResLayer2_mask]

#data from layers 1 and 2 separated 
Por1 = Porosity[ResLayer1_mask]
Por2 = Porosity[ResLayer2_mask]
Dens1 = Density[ResLayer1_mask]
Dens2 = Density[ResLayer2_mask]
Perm1 = Permeability[ResLayer1_mask]
Perm2 = Permeability[ResLayer2_mask]

#create porosity-permeability model for both layers using log-log
Por1_log   = log.(Por1[Perm1 .> 0.0])
Perm1_log  = log.(Perm1[Perm1 .> 0.0])
Por2_log   = log.(Por2[Perm2 .> 0.0])
Perm2_log  = log.(Perm2[Perm2 .> 0.0])

PorPerm1 = DataFrame(n=Por1_log, k=Perm1_log)
PorPerm2 = DataFrame(n=Por2_log, k=Perm2_log)
Formula2 = @formula(k ~ 1 + n)
LRPP1    = lm(Formula2, PorPerm1)
LRPP2    = lm(Formula2, PorPerm2)

#get regression line values for porosity-permeability plot
n_plot1 = range(minimum(Por1_log), maximum(Por1_log), length = 100)
PorPerm1_pred = predict(LRPP1, DataFrame(n = n_plot1), interval = :prediction, level = 0.95)
n_plot2 = range(minimum(Por2_log), maximum(Por2_log), length = 100)
PorPerm2_pred = predict(LRPP2, DataFrame(n = n_plot2), interval = :prediction, level = 0.95)

#extract prediction results 
predictionPP1 = Vector{Float64}(PorPerm1_pred.prediction)
lowerPP1      = Vector{Float64}(PorPerm1_pred.lower)
upperPP1      = Vector{Float64}(PorPerm1_pred.upper)
predictionPP2 = Vector{Float64}(PorPerm2_pred.prediction)
lowerPP2      = Vector{Float64}(PorPerm2_pred.lower)
upperPP2      = Vector{Float64}(PorPerm2_pred.upper)

#make plot
f = Figure()
ax1 = Axis(f[1, 1], 
title  = "Linear Regression 2 With Error Bounds", 
xlabel = "log(Porosity)",
ylabel = "log(Permeability)")

#plot for porosity/permeability
band!(ax1, n_plot1, lowerPP1, upperPP1, alpha = 0.3)
GLMakie.scatter!(ax1, Por1_log, Perm1_log, markersize=4)
lines!(ax1, n_plot1, predictionPP1, linewidth = 2)
lines!(ax1, n_plot1, lowerPP1, linewidth = 1.5, linestyle =:dash)
lines!(ax1, n_plot1, upperPP1, linewidth = 1.5, linestyle =:dash)

band!(ax1, n_plot2, lowerPP2, upperPP2, alpha = 0.3)
GLMakie.scatter!(ax1, Por2_log, Perm2_log, markersize=4)
lines!(ax1, n_plot2, predictionPP2, linewidth = 2)
lines!(ax1, n_plot2, lowerPP2, linewidth = 1.5, linestyle =:dash)
lines!(ax1, n_plot2, upperPP2, linewidth = 1.5, linestyle =:dash)

display(f)

I am aware of “for loops” but am not sure how I would go about implementing one in the code. Can anyone help?

2 
# First define a struct to hold your results here
struct Results
  # Exercise for reader
end
results = Results[]
for ilayer=1:2
  mask = LayerNum .== ilayer
  lx = XCoord[mask]
  ly = YCoord[mask]
  # etc.
  push!(results,Results(#=your regression data here=#))
end
# Plot or save results here
3

This would be a good place for a function. Something like this:

function compute_prediction(ResData)
    lx = ResData.X
    ly = ResData.Y
    por =  ResData.Porosity
    perm = ResData.Permeability
    Por_log  = log.(por[perm .> 0])
    Perm_log = log.(perm[perm .> 0])
    # ... keep going

    # return predictionPP, lowerPP, upperPP # can return as Tuple
    return (; predictionPP, lowerPP, upperPP) # but a NamedTuple labels the outputs nicely
end

# compute layers chosen manually
predictionPP1, lowerPP1, upperPP1 = compute_prediction(ResData[ResData.Layer .== 1, :]) # only processes where Layer == 1
predictionPP2, lowerPP2, upperPP2 = compute_prediction(ResData[ResData.Layer .== 2, :]) # only processes where Layer == 2

# or can compute a range of layers in a single call
results = map(layer -> compute_prediction(ResData[ResData.Layer .== layer, :]), 1:2) # compute layers in range 1:2
predictionPP1, lowerPP1, upperPP1 = results[1] # can extract them this way (note: if you don't start your layer range from 1, results[1] will not correspond to Layer==1)

(There’s a chance I’ve made minor typos here. I have not attempted to run this because obviously I didn’t finish writing everything.)