added julia models-incomplete

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
+# ╔═╡ 166031a9-d200-48f2-93f7-0e535d0a78b7
-begin #creating synthetic data
+begin
-	#successfully tested with n=100_000. It recovered the values quite well.
+	n=1000
-	n = 1000
+	enrollment_ratio = rand([0.2,0.5,0.7,1,1.1],n)
-
+	marketing = rand([0,1,2,1,5,7],n)
-	#####################Independent Data Draws ##############
+
-	#= Elapsed time 
+	x = hcat(enrollment_ratio,marketing)
-	How long the trial actually takes
+	k = size(x,2)
-	It should occur over a range of roughly [0,5], mostly clustered around 1
+	
-	=#
+	g=2 #number of categories
-	elapsed_time = rand(Gamma(8,0.125),n)
+	category = rand(1:g,n)
-
+
-	#= Enrollment success
+	M = reshape(1:(k*g),g,k) .* 0.1
-	How many of the initially planned group they have enrolled.
+	B = M .+ randn(g,k)
-	=#
+
-	enrollment_pctg = rand(Gamma(4,0.12),n)
+	mu = -1 .+ sum(x .* B[category,:],dims=2)
-
+	p = logistic.(mu)
-	#= Competing brands for compound
+	y = rand.(Bernoulli.(p))
-	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
+# ╔═╡ 2dc7ff65-08ed-4234-b997-50d999e54b32
-begin
+M
-	################# DEPENDENT VARIABLE CREATION ###############
-	#HyperParameters for conclusion_status
 
-	#draws of hyperparameters for impact of indication category
+# ╔═╡ 2654dac7-6ed8-4625-bb62-07ef771869b6
-	alpha_indication_category = [1,-0.4,0.7,2,0.1,-1]
 
-	#hyperparameters for impact of current state
-	alpha_current_state = [1,0.8,1.2,-1.2]
 
+# ╔═╡ b9fb5248-fff9-4b8e-8801-71a797e4c31e
 
-	#### 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)
+# ╔═╡ ec27451b-252f-406a-8d60-376f4dc79da1
 
-	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
+# ╔═╡ ef5d40d7-a9a8-40e4-8004-f40dbd58f837
-	p = logistic.(B)
-	y = rand.(Bernoulli.(p))
-end
 
-# ╔═╡ 322c7879-faf6-4f03-af2c-61e9c2cd76e3
-sum(y)
 
-# ╔═╡ 8b94ef35-4ec0-420a-9e2b-7aecf1497e97
+# ╔═╡ 0100a4ef-45d6-4ea3-b27c-93889bc09827
-histogram(p,bins = [-0.10,0.01,0.25,0.5,0.75,0.99,1.1])
 
-# ╔═╡ ea4df00b-72e2-4de6-afb1-f0b1d7e15e9c
-histogram(β_enrollment)
 
-# ╔═╡ 0696416b-72ba-46a0-ac54-0f8afbba554b
+# ╔═╡ 69e8107d-88d1-42d9-9a9b-93dbbf28bbe8
-@model function Logit(y,X,Z,enrollment_pctg)
+histogram(p,bins=20)
-	n,k_x = size(X)
-	k_z = size(Z,2)
 
-	α ~ MvNormal(zeros(10),I)
+# ╔═╡ 069c6761-ed63-4a16-aeaf-2fc4db029dfe
-	M = Z * α
+histogram(y)
-	#σ_β ~ Exponential(0.5)
-	#β_enrollment .~ Normal.(M,σ_β)
 
-	σ_everything_else ~ Exponential(0.5)
+# ╔═╡ be6c4cd5-a45d-42ea-845b-d7815c1e93f8
-	β_everything_else ~ MvNormal(zeros(k_x),σ_everything_else*I)
+@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
 
-	μ = X * β_everything_else .+ enrollment_pctg .* M
+	σ ~ filldist(Exponential(0.2),g)
 	
-	p = logistic.(μ)
+	β = Vector{Vector{T}}(undef,g)
-	y .~ Bernoulli.(p)
+	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
+# ╔═╡ 630d8bca-346d-49df-b62a-ebdd7530dd83
-testmodel = Logit(y,X,Z,enrollment_pctg);
+
+
+# ╔═╡ 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
+# ╔═╡ c518e08f-0a34-4622-bcf4-24e3f2235ed0
-chain = sample(testmodel,NUTS(),MCMCThreads(),2000,2)
 
-# ╔═╡ 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
+# ╠═166031a9-d200-48f2-93f7-0e535d0a78b7
-# ╠═75c5e5a8-0fe5-4739-99a3-4cc17487dc91
+# ╠═2dc7ff65-08ed-4234-b997-50d999e54b32
-# ╠═322c7879-faf6-4f03-af2c-61e9c2cd76e3
+# ╠═2654dac7-6ed8-4625-bb62-07ef771869b6
-# ╠═8b94ef35-4ec0-420a-9e2b-7aecf1497e97
+# ╠═b9fb5248-fff9-4b8e-8801-71a797e4c31e
-# ╠═ea4df00b-72e2-4de6-afb1-f0b1d7e15e9c
+# ╠═ec27451b-252f-406a-8d60-376f4dc79da1
-# ╠═0696416b-72ba-46a0-ac54-0f8afbba554b
+# ╠═ef5d40d7-a9a8-40e4-8004-f40dbd58f837
-# ╠═3f4aa70d-ffce-44fe-8520-831a5a94fe47
+# ╠═0100a4ef-45d6-4ea3-b27c-93889bc09827
-# ╠═989b69de-e187-48f2-8454-8d86b232d895
+# ╠═69e8107d-88d1-42d9-9a9b-93dbbf28bbe8
-# ╠═8da7994b-ac63-4583-a360-d49cc8b04d35
+# ╠═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