diff --git a/Julia-turing models/PlutoModel.jl b/Julia-turing models/PlutoModel.jl index ebdc478..5b56c12 100644 --- a/Julia-turing models/PlutoModel.jl +++ b/Julia-turing models/PlutoModel.jl @@ -18,123 +18,103 @@ md""" So it seems like I am having issues where I can't get the herarcheal parameters to work. """ -# ╔═╡ bd9391b7-5240-46ca-be26-357292ffaddb -begin #creating synthetic data - #successfully tested with n=100_000. It recovered the values quite well. - n = 1000 - - #####################Independent Data Draws ############## - #= Elapsed time - How long the trial actually takes - It should occur over a range of roughly [0,5], mostly clustered around 1 - =# - elapsed_time = rand(Gamma(8,0.125),n) - - #= Enrollment success - How many of the initially planned group they have enrolled. - =# - enrollment_pctg = rand(Gamma(4,0.12),n) - - #= Competing brands for compound - The number of other brands selling that compound. - =# - compound_count = rand(Poisson(2),n) - - #= Indication Category - Ignore for now - =# - indication_category_index = rand(1:6,n) - indication_category_dummies = zeros(n,6) - for (row,col) in enumerate(indication_category_index) - indication_category_dummies[row,col] = 1 - end - #= Indication Population rate - measured in per 10_000, across 6 national develoment measures - =# - population_rate_raw = rand(MvNormal([50, 60,100, 70,65,40],5*I),n) - population_rate = population_rate_raw[indication_category_index] - - #= Current State - One of 3 options - - Recruiting - - Not yet recruiting - - Active, not recruiting - - Suspended - =# - current_state_prob = Categorical([0.5,0.05,0.4,0.05]) - current_state_index = rand(current_state_prob,n) - current_state_dummies = zeros(n,4) - for (row,col) in enumerate(current_state_index) - current_state_dummies[row,col] = 1 - end - Z = hcat(indication_category_dummies,current_state_dummies) - X = hcat(enrollment_pctg,compound_count,population_rate) -end - -# ╔═╡ 75c5e5a8-0fe5-4739-99a3-4cc17487dc91 +# ╔═╡ 166031a9-d200-48f2-93f7-0e535d0a78b7 begin - ################# DEPENDENT VARIABLE CREATION ############### - #HyperParameters for conclusion_status + n=1000 + enrollment_ratio = rand([0.2,0.5,0.7,1,1.1],n) + marketing = rand([0,1,2,1,5,7],n) + + x = hcat(enrollment_ratio,marketing) + k = size(x,2) + + g=2 #number of categories + category = rand(1:g,n) + + M = reshape(1:(k*g),g,k) .* 0.1 + B = M .+ randn(g,k) - #draws of hyperparameters for impact of indication category - alpha_indication_category = [1,-0.4,0.7,2,0.1,-1] + mu = -1 .+ sum(x .* B[category,:],dims=2) + p = logistic.(mu) + y = rand.(Bernoulli.(p)) +end - #hyperparameters for impact of current state - alpha_current_state = [1,0.8,1.2,-1.2] +# ╔═╡ 2dc7ff65-08ed-4234-b997-50d999e54b32 +M +# ╔═╡ 2654dac7-6ed8-4625-bb62-07ef771869b6 - #### Parameters for conclusion_status - #get β for enrollment - M_enrollment = Z * vcat(alpha_indication_category,alpha_current_state) - S_enrollment = 0.7 - β_enrollment = M_enrollment .+ S_enrollment .* randn(n) +# ╔═╡ b9fb5248-fff9-4b8e-8801-71a797e4c31e - k_x = size(X,2) - β_everything_else = rand(MvNormal([0.5,0.4,0.3],I),n) - B = sum(X .* β_everything_else',dims=2) .+ β_enrollment .* enrollment_pctg - p = logistic.(B) - y = rand.(Bernoulli.(p)) -end +# ╔═╡ ec27451b-252f-406a-8d60-376f4dc79da1 -# ╔═╡ 322c7879-faf6-4f03-af2c-61e9c2cd76e3 -sum(y) -# ╔═╡ 8b94ef35-4ec0-420a-9e2b-7aecf1497e97 -histogram(p,bins = [-0.10,0.01,0.25,0.5,0.75,0.99,1.1]) +# ╔═╡ ef5d40d7-a9a8-40e4-8004-f40dbd58f837 -# ╔═╡ ea4df00b-72e2-4de6-afb1-f0b1d7e15e9c -histogram(β_enrollment) -# ╔═╡ 0696416b-72ba-46a0-ac54-0f8afbba554b -@model function Logit(y,X,Z,enrollment_pctg) - n,k_x = size(X) - k_z = size(Z,2) +# ╔═╡ 0100a4ef-45d6-4ea3-b27c-93889bc09827 - α ~ MvNormal(zeros(10),I) - M = Z * α - #σ_β ~ Exponential(0.5) - #β_enrollment .~ Normal.(M,σ_β) - σ_everything_else ~ Exponential(0.5) - β_everything_else ~ MvNormal(zeros(k_x),σ_everything_else*I) +# ╔═╡ 69e8107d-88d1-42d9-9a9b-93dbbf28bbe8 +histogram(p,bins=20) - μ = X * β_everything_else .+ enrollment_pctg .* M +# ╔═╡ 069c6761-ed63-4a16-aeaf-2fc4db029dfe +histogram(y) - p = logistic.(μ) - y .~ Bernoulli.(p) +# ╔═╡ be6c4cd5-a45d-42ea-845b-d7815c1e93f8 +@model function hielogit(y, cat, x,g, ::Type{T} = Float64) where {T} + n,k = size(x) + k = k+1 + x1 = ones(n) + x = hcat(x1,x) + μ ~ MvNormal(k,1) #hyperparam + + σ ~ filldist(Exponential(0.2),g) + β = Vector{Vector{T}}(undef,g) + for i in 1:g + β[i] ~ MvNormal(μ,σ[i]) + end + + for i in 1:n + p = sum(x[1,:] .* β[cat[i]]) + y[i] ~ BernoulliLogit(p) + end end -# ╔═╡ 3f4aa70d-ffce-44fe-8520-831a5a94fe47 -testmodel = Logit(y,X,Z,enrollment_pctg); +# ╔═╡ 630d8bca-346d-49df-b62a-ebdd7530dd83 + + +# ╔═╡ e5c05073-7686-4053-8593-3953d8f18f16 + + +# ╔═╡ 153ae545-5baf-4c53-a76c-f22a72994096 + + +# ╔═╡ f113c911-64ab-4998-8a1c-7b17af4b2e9a +hl = hielogit(y,category,x,g) + +# ╔═╡ 37f3dadf-d07b-49c7-82fb-5ce49edc4476 +chains_of_jacob_marley = sample(hl, NUTS(),2000) + +# ╔═╡ 887d70f4-3592-4cd9-b290-66e07a4821f9 +plot(chains_of_jacob_marley) + +# ╔═╡ 09bf1ae5-3b22-4c40-a99b-7eab5aafd76d +chains_of_jacob_marley.name_map + +# ╔═╡ a4ba8cc8-e187-4841-ae2e-dbc94b1bfc76 +mean(chains_of_jacob_marley.value.data[:,1:11,1],dims=1) + +# ╔═╡ 25bf4f46-0051-4d06-8582-38dd4ac3ecdf + + +# ╔═╡ be600b5b-745b-4750-b534-a5a40461344a +chains_of_jacob_marley.value.axes[2].val[1:11] -# ╔═╡ 989b69de-e187-48f2-8454-8d86b232d895 -chain = sample(testmodel,NUTS(),MCMCThreads(),2000,2) +# ╔═╡ c518e08f-0a34-4622-bcf4-24e3f2235ed0 -# ╔═╡ 8da7994b-ac63-4583-a360-d49cc8b04d35 -plot(chain) # ╔═╡ 00000000-0000-0000-0000-000000000001 PLUTO_PROJECT_TOML_CONTENTS = """ @@ -161,7 +141,7 @@ PLUTO_MANIFEST_TOML_CONTENTS = """ julia_version = "1.8.5" manifest_format = "2.0" -project_hash = "796d12ee251d971e0e0bf81a91d64b9c58913e1d" +project_hash = "2a5d2f6d90e4bc9a5bca3804cd573e1dbc6b5b6e" [[deps.AbstractFFTs]] deps = ["ChainRulesCore", "LinearAlgebra"] @@ -1794,14 +1774,26 @@ version = "1.4.1+0" # ╠═54680695-4d31-4dd9-8b3c-c0af5153abef # ╠═735ae519-019e-48c6-8b73-7982fa60501a # ╠═f6fc4578-43c5-4fb7-a24c-2567aa5b33c1 -# ╠═bd9391b7-5240-46ca-be26-357292ffaddb -# ╠═75c5e5a8-0fe5-4739-99a3-4cc17487dc91 -# ╠═322c7879-faf6-4f03-af2c-61e9c2cd76e3 -# ╠═8b94ef35-4ec0-420a-9e2b-7aecf1497e97 -# ╠═ea4df00b-72e2-4de6-afb1-f0b1d7e15e9c -# ╠═0696416b-72ba-46a0-ac54-0f8afbba554b -# ╠═3f4aa70d-ffce-44fe-8520-831a5a94fe47 -# ╠═989b69de-e187-48f2-8454-8d86b232d895 -# ╠═8da7994b-ac63-4583-a360-d49cc8b04d35 +# ╠═166031a9-d200-48f2-93f7-0e535d0a78b7 +# ╠═2dc7ff65-08ed-4234-b997-50d999e54b32 +# ╠═2654dac7-6ed8-4625-bb62-07ef771869b6 +# ╠═b9fb5248-fff9-4b8e-8801-71a797e4c31e +# ╠═ec27451b-252f-406a-8d60-376f4dc79da1 +# ╠═ef5d40d7-a9a8-40e4-8004-f40dbd58f837 +# ╠═0100a4ef-45d6-4ea3-b27c-93889bc09827 +# ╠═69e8107d-88d1-42d9-9a9b-93dbbf28bbe8 +# ╠═069c6761-ed63-4a16-aeaf-2fc4db029dfe +# ╠═be6c4cd5-a45d-42ea-845b-d7815c1e93f8 +# ╠═630d8bca-346d-49df-b62a-ebdd7530dd83 +# ╠═e5c05073-7686-4053-8593-3953d8f18f16 +# ╠═153ae545-5baf-4c53-a76c-f22a72994096 +# ╠═f113c911-64ab-4998-8a1c-7b17af4b2e9a +# ╠═37f3dadf-d07b-49c7-82fb-5ce49edc4476 +# ╠═887d70f4-3592-4cd9-b290-66e07a4821f9 +# ╠═09bf1ae5-3b22-4c40-a99b-7eab5aafd76d +# ╠═a4ba8cc8-e187-4841-ae2e-dbc94b1bfc76 +# ╠═25bf4f46-0051-4d06-8582-38dd4ac3ecdf +# ╠═be600b5b-745b-4750-b534-a5a40461344a +# ╠═c518e08f-0a34-4622-bcf4-24e3f2235ed0 # ╟─00000000-0000-0000-0000-000000000001 # ╟─00000000-0000-0000-0000-000000000002