Create a dataset used in a marginal structural model will all values of mixture set to a specific value and covariates (if any) set to reference values.

Confidence bounds for effect measures and linear predictions at joint exposure values

using Qgcomp, DataFrames, StatsModels

x1 = rand(100, 3)
x = rand(100, 3)
z = rand(100, 3)
xq, _ = Qgcomp.get_xq(x, 4)
y = randn(100) + xq * [.1, 0.05, 0]
data = DataFrame(hcat(y,x,z), [:y, :x1, :x2, :x3, :z1, :z2, :z3])
form = @formula(y~x1+x2+x3+z1+z2+z3)
form_noint = @formula(y~ -1+x1+x2+x3+z1+z2+z3)
expnms = [:x1, :x2, :x3]

m = qgcomp_glm_ee(form, data, expnms, 4, Normal())
Qgcomp.printbounds(bounds(m, 0:0.1:3, 0.8))

bsplineIntknts(0,5, 8)

bsplineZeroBasis([0.0,0.3,3.0],0.1, 2.0)

bsplineknots(0,5, 8) extdist = determines the distance between external knots. Set to 0 for classic "repeated" tail knots. if positive, a constant gap between exterior knots starting at max value + epsilonilon if negative, a constant multiplier of the final gap size between interior knots

using Qgcomp, DataFrames, Distributions, StatsModels

n = 200
dat = DataFrame(y=Int.(rand(Bernoulli(0.25), n)), x1=rand(n), x2=rand(n), z=rand(n))

# 3 internal knots
deg = 2
knts = bsplineknots(dat.x1, 5, deg; ptype = "uniform", extdist = 1)
nkn = length(knts)
nkn - deg 
size(bs(dat.x1,knts), 2)


kform = term(:y) ~ term(1) + bs(:x1,knts) + term(:x2)
ft = qgcomp_glm_ee(kform, dat,[:x1, :x2], nothing, Normal())
ft = qgcomp_glm_boot(kform, dat,[:x1, :x2], nothing, Normal())


using GLM
glm(kform, dat, Binomial())

using Qgcomp
using LSurvival, DataFrames, Random
rng = MersenneTwister()
# expected effect size in qgcomp for X1, X2 ≈ (truebeta[1] + truebeta[2])/4 = 0.5
truebeta = [4.0, -2.0, 1.0, -1.0, 1.0]
approxpsi = (truebeta[1] + truebeta[2])/4
X, t, d = LSurvival.dgm_phmodel(300; λ=1.25,β=truebeta)
survdata = hcat(DataFrame(X, [:x1, :x2, :z1, :z2, :z3]), DataFrame(hcat(t,d),[:t,:d]))

rng = Xoshiro(1232)

# conditional MSM with fast estimator
qgcomp_cox_noboot(rng, @formula(Surv(t, d)~x1+x2+z1+z2+z3), survdata, ["x1", "x2"], 4)

# conditional MSM with traditional g-computation estimator (conditional on covariates - should look a lot like the "noboot" version)
qgcomp_cox_boot(rng, @formula(Surv(t, d)~x1+x2+z1+z2+z3), survdata, ["x1", "x2"], 4, msmformula=@formula(Surv(t, d)~mixture+z1+z2+z3))

# population MSM with traditional g-computation estimator (necessary, in this case)
qgcomp_cox_boot(rng, @formula(Surv(t, d)~x1+x2+z1+z2+z3), survdata, ["x1", "x2"], 4, msmformula=@formula(Surv(t, d)~mixture+z1+z2+z3))

# non-linear MSM with traditional g-computation estimator (necessary, in this case)
qgcomp_cox_boot(rng, @formula(Surv(t, d)~x1+x2+x1*x2+z1+z2+z3), survdata, ["x1", "x2"], 4, msmformula=@formula(Surv(t, d)~mixture+mixture^2))

using Qgcomp
using LSurvival, DataFrames, Random
id, int, out, data = LSurvival.dgm(MersenneTwister(1212), 100, 20);

data[:, 1] = round.(data[:, 1], digits = 3);
d, X = data[:, 4], data[:, 1:3];
wt = ones(length(d)) # random weights just to demonstrate usage
tab = ( in = int, out = out, d=d, x=X[:,1], z1=X[:,2], z2=X[:,3]) ;
df = DataFrame(tab)

coxph(@formula(Surv(in, out, d)~x+z1+z2), tab, ties = "efron", wts = wt) |> display
m = qgcomp_cox_noboot(@formula(Surv(in, out, d)~x+z1+z2), df, ["z1", "z2"], 4)

using Qgcomp, DataFrames, StatsModels, StatsBase

x1 = rand(100, 3)
x = rand(100, 3)
z = rand(100, 3)
xq, _ = Qgcomp.get_xq(x, 4)
y = randn(100) + xq * [.1, 0.05, 0]+ (xq .* xq) * [-.1, 0.05, 0]
ybin = Int.(y .> median(y))
data = DataFrame(hcat(ybin,y,x,z), [:ybin, :y, :x1, :x2, :x3, :z1, :z2, :z3])
form = @formula(y~x1+x2+x3+x1^2+x2^2+x3^2+z1+z2+z3)
formbin = @formula(ybin~x1+x2+x3+x1^2+x2^2+x3^2+z1+z2+z3)
expnms = [:x1, :x2, :x3]

# note the top fit is incorrect
m0 = qgcomp_glm_noboot(form, data, expnms, 4, Normal()) 

# three ways to specify non-linear fits
m = qgcomp_glm_boot(form, data, expnms, 4, Normal(), B=2000, msmformula=@formula(y~mixture+mixture^2)) 
mb = qgcomp_glm_boot(form, data, expnms, 4, Normal(), B=2000, degree=2) 
m2 = qgcomp_glm_ee(form, data, expnms, 4, Normal(), degree=2) 
isfitted(m)
fitted(m)
aic(m)
aicc(m)
bic(m)
loglikelihood(m)

# binary outcome
# note the top fit is incorrect
m0 = qgcomp_glm_noboot(formbin, data, expnms, 4, Bernoulli()) 

# three ways to specify non-linear fits
m = qgcomp_glm_boot(formbin, data, expnms, 4, Binomial(), B=2000, msmformula=@formula(y~mixture+mixture^2)) 
mb = qgcomp_glm_boot(formbin, data, expnms, 4, Binomial(), B=2000, degree=2) 
m2 = qgcomp_glm_ee(formbin, data, expnms2, 4, Binomial(), degree=2) 

binary

dat = DataFrame(y=Int.(rand(Bernoulli(0.25), 50)), x1=rand(50), x2=rand(50), z=rand(50))

# Marginal mixture OR (no confounders)
ft1 = qgcomp_glm_noboot(@formula(y ~ x1 + x2), dat,["x1", "x2"], nothing, Binomial())
ft2 = qgcomp_glm_ee(@formula(y ~ x1 + x2), dat,["x1", "x2"], nothing, Binomial())
ft3 = qgcomp_glm_ee(@formula(y ~ x1 + x2), dat,["x1", "x2"], nothing, Binomial(), rr=true)

# Conditional mixture OR
qgcomp_glm_noboot(@formula(y ~ z + x1 + x2), dat,["x1", "x2"], 4, Binomial())
# Marginal mixture OR
qgcomp_glm_ee(@formula(y ~ z + x1 + x2), dat,["x1", "x2"], 4, Binomial())

using Qgcomp, DataFrames, StatsModels

x1 = rand(100, 3)
x = rand(100, 3)
z = rand(100, 3)
xq, _ = Qgcomp.get_xq(x, 4)
y = randn(100) + xq * [.1, 0.05, 0]
data = DataFrame(hcat(y,x,z), [:y, :x1, :x2, :x3, :z1, :z2, :z3])
form = @formula(y~x1+x2+x3+z1+z2+z3)
form_noint = @formula(y~-1+x1+x2+x3+z1+z2+z3)
expnms = [:x1, :x2, :x3]

m = qgcomp_glm_noboot(form, data, expnms, 4, Normal())
m = qgcomp_glm_noboot(form_noint, data, expnms, 4, Normal())
fitted(m)
aic(m)
aicc(m)
bic(m)
loglikelihood(m)

x = rand(100, 3)
q=4
nexp = size(x,2)

xq, breaks = quantize(x, q)

Create a vector of restricted spline knots (based on rcspline.eval from R Hmisc package, normalized version) usage: rsplineknots(x,nk) x = vector of numbers nk = number of (interior) knots output: a nk-length vector of (interior) knots

```julia using Qgcomp, DataFrames, Distributions, StatsModels dat = DataFrame(y=Int.(rand(Bernoulli(0.25), 50)), x1=rand(50), x2=rand(50), z=rand(50))

knts = rsplineknots(dat.x1, 5)

kform = term(:y)~ term(1) + rcs(:x1,knts) + term(:x2) kform2 = term(:y)~ term(1) + rqs(:x1,knts) + term(:x2)

ft = qgcompglmee(kform, dat,["x1", "x2"], nothing, Binomial()) ft = qgcompglmee(kform2, dat,["x1", "x2"], nothing, Binomial())

```

Covariance between linear combinations of terms

e.g. for a linear regression β0 + β1x + β2z + β3w + β4r

with coefficient covariance matrix V,

the variance for ψ = β1 + β2 (e.g. the simultaneous effect of increasing x and z by one unit) is calculated as vccomb(V, [0,1,1,0,0], [0,1,1,0,0])

The covariance term COV(ψ, β0) is calculated as vccomb(V, [0,1,1,0,0], [1,0,0,0,0])

This function is useful for qgcomp methods when ψ is just a sum of β coefficients, because it allows straightforward calculation of the full covariance matrix. The underlying calculations are very simple, but the function provides

One version of the function (string vector arguments) is very specifically tuned toward a covariance matrix for qgcompglmboot, given the covariance matrix column labels and expnms. It assumes a univariate ψ parameter.

regression curve plot for msm

using Qgcomp
using Random, GLM, DataFrames, LSurvival
using Plots

#using LinearAlgebra, StatsBase
pointwise = Qgcomp.pointwise
#loesspred = Qgcomp.loesspred
smoothpred = Qgcomp.loesspred

x1 = rand(100, 3)
x = rand(100, 3)
z = rand(100, 3)
xq, _ = Qgcomp.get_xq(x, 4)
y = randn(100) + xq .*  xq * [0.1, 0.05, 0]
lindata = DataFrame(hcat(y, x, z), [:y, :x1, :x2, :x3, :z1, :z2, :z3])

m = qgcomp_glm_boot(@formula(y~x1+x1^2+x2+x2^2+x3+x3^2+z1+z2+z3), lindata, ["x1", "x2", "x3", "z1"], 6, Normal(), msmformula=@formula(y~mixture+mixture^2))
#x,y,ylower,yupper = pointwise(m.msm)

responseplot(m)

# plots based on these bounds, except the smooth function
bnds = bounds(m)
mw = bnds[:model]
pw = bnds[:pointwise]

Plotting weights

using Qgcomp
using Random, GLM, DataFrames, LSurvival
using Plots

x1 = rand(100, 3)
x = rand(100, 3)
z = rand(100, 3)
xq, _ = Qgcomp.get_xq(x, 4)
y = randn(100) + xq * [0.1, 0.05, 0]
lindata = DataFrame(hcat(y, x, z), [:y, :x1, :x2, :x3, :z1, :z2, :z3])

mint = qgcomp_glm_noboot(@formula(y~x1+x2+x3+z1+z2+z3), lindata, ["x1", "x2", "x3", "z1", "z2", "z3"], 4, Normal())
weightplot(mint)

n=300
x = rand(n, 3)
z = rand(n, 3)
xq, _ = Qgcomp.get_xq(x, 4)
y = randn(n) + xq * [.1, .05, 0]
data = DataFrame(hcat(y,x,z), [:y, :x1, :x2, :x3, :z1, :z2, :z3])
form = @formula(y~x1+x2+x3+x3+z1+z2+z3)
form2 = @formula(y~x1+x2+x3+x1^2+x2^2+x3^2+x1*z1+z2+z3)
msmformula = @formula(y~mixture+mixture^2)
expnms = [:x1, :x2, :x3]
q = 4
m = QGcomp_glm(form, data, expnms, 4, Normal());
fit!(m)
m2 = QGcomp_ee(form, data, expnms, 4, Normal());
StatsBase.fit!(m2)
m = QGcomp_ee(form, data, expnms, q, Normal())
qgcomp_glm_noboot(form, data, expnms, 4, Normal())
ft = qgcomp_glm_ee(form2, data, expnms, q, Normal(),degree=2)



keys = Dict([
    :msmformula => msmformula,
])
contrasts = Dict{Symbol,Any}()
rr = false

using Qgcomp, DataFrames, StatsModels

x = rand(100, 3)
z = rand(100, 3)
xq, _ = Qgcomp.get_xq(x, 4)
y = randn(100) + xq * [.1, .05, 0]
data = DataFrame(hcat(y,x,z), [:y, :x1, :x2, :x3, :z1, :z2, :z3])
formula = @formula(y~x1+x2+x3+z1+z2+z3)
expnms = ["x"*string(i) for i in 1:3]
m = Qgcomp.QGcomp_glm(formula, data, expnms, 4, Normal())
m 
fit!(m)
m