diff --git a/docs/src/fig/exponential.svg b/docs/src/fig/exponential.svg deleted file mode 100644 index e56ced7..0000000 --- a/docs/src/fig/exponential.svg +++ /dev/null @@ -1,38 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/docs/src/fig/exponential_haz.svg b/docs/src/fig/exponential_haz.svg new file mode 100644 index 0000000..7b5de66 --- /dev/null +++ b/docs/src/fig/exponential_haz.svg @@ -0,0 +1,38 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/fig/gamma.svg b/docs/src/fig/gamma.svg deleted file mode 100644 index dcb20e5..0000000 --- a/docs/src/fig/gamma.svg +++ /dev/null @@ -1,38 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/docs/src/fig/gamma_haz.svg b/docs/src/fig/gamma_haz.svg new file mode 100644 index 0000000..cd38e8f --- /dev/null +++ b/docs/src/fig/gamma_haz.svg @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/fig/gengamma_haz.svg b/docs/src/fig/gengamma_haz.svg new file mode 100644 index 0000000..9cb44a1 --- /dev/null +++ b/docs/src/fig/gengamma_haz.svg @@ -0,0 +1,38 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/fig/gengamma_pdf.svg b/docs/src/fig/gengamma_pdf.svg deleted file mode 100644 index 44b4edd..0000000 --- a/docs/src/fig/gengamma_pdf.svg +++ /dev/null @@ -1,39 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/docs/src/fig/lognormal_haz.svg b/docs/src/fig/lognormal_haz.svg new file mode 100644 index 0000000..27bc163 --- /dev/null +++ b/docs/src/fig/lognormal_haz.svg @@ -0,0 +1,39 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/fig/lognormal_pdf.svg b/docs/src/fig/lognormal_pdf.svg deleted file mode 100644 index 0f239d0..0000000 --- a/docs/src/fig/lognormal_pdf.svg +++ /dev/null @@ -1,37 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/docs/src/fig/weibull_haz.svg b/docs/src/fig/weibull_haz.svg new file mode 100644 index 0000000..9a22359 --- /dev/null +++ b/docs/src/fig/weibull_haz.svg @@ -0,0 +1,37 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/fig/survdist_pdf.svg b/docs/src/fig/weibull_pdf.svg similarity index 85% rename from docs/src/fig/survdist_pdf.svg rename to docs/src/fig/weibull_pdf.svg index fcee2f4..9fb2a22 100644 --- a/docs/src/fig/survdist_pdf.svg +++ b/docs/src/fig/weibull_pdf.svg @@ -1,39 +1,39 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/fig/weibull_risk.svg b/docs/src/fig/weibull_risk.svg new file mode 100644 index 0000000..ea42f1b --- /dev/null +++ b/docs/src/fig/weibull_risk.svg @@ -0,0 +1,38 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/fig/survdist_surv.svg b/docs/src/fig/weibull_surv.svg similarity index 85% rename from docs/src/fig/survdist_surv.svg rename to docs/src/fig/weibull_surv.svg index 333ef48..466e8ea 100644 --- a/docs/src/fig/survdist_surv.svg +++ b/docs/src/fig/weibull_surv.svg @@ -1,38 +1,38 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/parametric.md b/docs/src/parametric.md index e12340c..44480f8 100644 --- a/docs/src/parametric.md +++ b/docs/src/parametric.md @@ -89,36 +89,53 @@ using Plots aftdist(weibullfit, label="X=0, Z1=0, Z2=0", title="Weibull distribution") # for specific covariate levels, include a 1.0 for the intercept aftdist!(weibullfit, covlevels=[1.0, 1.0, 1.0, 1.0], color="red", label="X=1, Z1=1, Z2=1", npoints=300) -savefig("survdist_pdf.svg") +savefig("weibull_pdf.svg") ``` -![Weibull-density](fig/survdist_pdf.svg) +![Weibull-density](fig/weibull_pdf.svg) ## Visualizing the distributions, survival distribution ```julia -using Plots - aftdist(weibullfit, type="surv", label="X=0, Z1=0, Z2=0", title="Weibull distribution") # for specific covariate levels, include a 1.0 for the intercept aftdist!(weibullfit, type="surv", covlevels=[1.0, 1.0, 1.0, 1.0], color="red", label="X=1, Z1=1, Z2=1") -savefig("survdist_surv.svg") +savefig("weibull_surv.svg") +``` +![Weibull-survival](fig/weibull_surv.svg) + +## Visualizing the distributions, hazard function +```julia +aftdist(weibullfit, type="haz", label="X=0, Z1=0, Z2=0", title="Weibull distribution") +# for specific covariate levels, include a 1.0 for the intercept +aftdist!(weibullfit, type="haz", covlevels=[1.0, 1.0, 1.0, 1.0], color="red", label="X=1, Z1=1, Z2=1", npoints=300) +savefig("weibull_haz.svg") ``` -![Weibull-survival](fig/survdist_surv.svg) +![Weibull-haz](fig/weibull_haz.svg) + +## Visualizing the distributions, risk function +```julia +aftdist(weibullfit, type="risk", label="X=0, Z1=0, Z2=0", title="Weibull distribution") +# for specific covariate levels, include a 1.0 for the intercept +aftdist!(weibullfit, type="risk", covlevels=[1.0, 1.0, 1.0, 1.0], color="red", label="X=1, Z1=1, Z2=1", npoints=300) +savefig("weibull_risk.svg") +``` +![Weibull-risk](fig/weibull_risk.svg) + ## Other distributions ### Exponential -Note that the exponential fit looks a lot like the Weibull fit. The log-scale parameter in the Weibull fit is close to zero. When it is exactly zero, then the Weibull distribution and exponential distribution are identical. +Note that the exponential fit assumes a constant hazard. ```julia expfit = survreg(@formula(Surv(in, out, d)~x+z1+z2), tab, wts=tab.wts, dist=LSurvival.Exponential()) -aftdist(expfit, type="surv", label="X=0, Z1=0, Z2=0", title="Exponential distribution") +aftdist(expfit, type="haz", label="X=0, Z1=0, Z2=0", title="Exponential distribution") # for specific covariate levels, include a 1.0 for the intercept -aftdist!(expfit, type="surv", covlevels=[1.0, 1.0, 1.0, 1.0], color="red", label="X=1, Z1=1, Z2=1") -savefig("exponential.svg") +aftdist!(expfit, type="haz", covlevels=[1.0, 1.0, 1.0, 1.0], color="red", label="X=1, Z1=1, Z2=1") +savefig("exponential_haz.svg") ``` -![exponential](fig/exponential.svg) +![exponential](fig/exponential_haz.svg) Output: @@ -149,11 +166,11 @@ Solver iterations: 14 # Here are results from a simpler model dat1 = (time = [1, 1, 6, 6, 8, 9], status = [1, 0, 1, 1, 0, 1], x = [1, 1, 1, 0, 0, 0]) lognormalfit = survreg(@formula(Surv( time, status)~x), dat1, dist=LSurvival.Lognormal()) -aftdist(lognormalfit, label="X=0", title="Log-normal distribution") -aftdist!(lognormalfit, covlevels=[1.0], color="red", label="X=1") -savefig("lognormal_pdf.svg") +aftdist(lognormalfit, type="haz", label="X=0", title="Log-normal distribution") +aftdist!(lognormalfit, type="haz", covlevels=[1.0], color="red", label="X=1") +savefig("lognormal_haz.svg") ``` -![log-normal](fig/lognormal_pdf.svg) +![log-normal](fig/lognormal_haz.svg) ```output Maximum likelihood estimates (alpha=0.05): @@ -175,12 +192,12 @@ Solver iterations: 9 ### Gamma ```julia gammafit = survreg(@formula(Surv(in, out, d)~x+z1+z2), tab, wts=tab.wts, dist=LSurvival.Gamma()) -aftdist(gammafit, type="surv", label="X=0, Z1=0, Z2=0", title="Gamma distribution") +aftdist(gammafit, type="haz", label="X=0, Z1=0, Z2=0", title="Gamma distribution") # for specific covariate levels, include a 1.0 for the intercept -aftdist!(gammafit, type="surv", covlevels=[1.0, 1.0, 1.0, 1.0], color="red", label="X=1, Z1=1, Z2=1") -savefig("gamma.svg") +aftdist!(gammafit, type="haz", covlevels=[1.0, 1.0, 1.0, 1.0], color="red", label="X=1, Z1=1, Z2=1") +savefig("gamma_haz.svg") ``` -![gamma](fig/gamma.svg) +![gamma](fig/gamma_haz.svg) @@ -221,11 +238,11 @@ ggammafit2 = survreg(@formula(Surv(in, out, d)~x+z1+z2), tab, wts=tab.wts, dist= ) ggammafit2 = survreg(@formula(Surv(t, d)~x), wtab, dist=LSurvival.GGamma()) -aftdist(ggammafit2, label="X=0", title="Generalized gamma distribution") -aftdist!(ggammafit2, covlevels=[1.0], color="red", label="X=1") -savefig("gengamma_pdf.svg") +aftdist(ggammafit2, type="haz", label="X=0", title="Generalized gamma distribution") +aftdist!(ggammafit2, type="haz", covlevels=[1.0], color="red", label="X=1") +savefig("gengamma_haz.svg") ``` -![generalized-gamma](fig/gengamma_pdf.svg) +![generalized-gamma](fig/gengamma_haz.svg) ```output ┌ Warning: Optimizer reports model did not converge. Gradient: [-0.30185255746883194, -92.76305273476676, -12.51874622578631] diff --git a/src/plot_recipes.jl b/src/plot_recipes.jl index 003db1b..dac3c8a 100644 --- a/src/plot_recipes.jl +++ b/src/plot_recipes.jl @@ -368,9 +368,22 @@ dat2 = ( ) fte = survreg(@formula(Surv(enter, exit, status)~x), dat2) +# density function aftdist(fte, label="X=0") aftdist!(fte, covlevels=[1.0, 2.0], color="red", label="X=1") +# Survival function +aftdist(fte, type="surv", label="X=0") +aftdist!(fte, type="surv", covlevels=[1.0, 2.0], color="red", label="X=1") + +# hazard function +aftdist(fte, type="haz", label="X=0") +aftdist!(fte, type="haz", covlevels=[1.0, 2.0], color="red", label="X=1") + +# Cumulative incidence/risk function +aftdist(fte, type="risk", label="X=0") +aftdist!(fte, type="risk", covlevels=[1.0, 2.0], color="red", label="X=1") + ``` """ @recipe function f(h::AFTdist; type = "pdf", covlevels = nothing, npoints = 100) @@ -389,14 +402,22 @@ aftdist!(fte, covlevels=[1.0, 2.0], color="red", label="X=1") end plotdist = name(typeof(dist))(sum(coefs .* covlevels), ft.P._S...) times = range(timeminmax[1], timeminmax[2], npoints) - if type == "pdf" + if type[1:3] == "pdf" dist = [exp(lpdf(plotdist, t)) for t in times] ylab --> "Density" - elseif type == "surv" + elseif type[1:3] == "sur" dist = [exp(lsurv(plotdist, t)) for t in times] ylab --> "Survival" + elseif type[1:3] == "ris" + dist = 1 .- [exp(lsurv(plotdist, t)) for t in times] + ylab --> "Risk" + elseif type[1:3] == "haz" + S = [exp(lsurv(plotdist, t)) for t in times] + f = [exp(lpdf(plotdist, t)) for t in times] + dist = f ./ S + ylab --> "Hazard" else - throw("Type must either be 'surv' or 'pdf'") + throw("Type must either be 'surv', 'risk', 'haz' or 'pdf'") end xlab --> "Time" grid --> false