Skip to content

Commit

Permalink
minor change
Browse files Browse the repository at this point in the history
  • Loading branch information
@IgorDouven
IgorDouven committed Dec 26, 2020
1 parent 182ab66 commit 8a6edfb
Show file tree
Hide file tree
Showing 2 changed files with 389 additions and 144 deletions.
354 changes: 354 additions & 0 deletions .ipynb_checkpoints/Ch3-checkpoint.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,354 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using DataFrames\n",
"using MixedModels\n",
"using StatsBase\n",
"using Gadfly\n",
"using CSV\n",
"using CategoricalArrays\n",
"using Colors"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"function gen_brew_colors(n)\n",
" cs = distinguishable_colors(n, \n",
" [colorant\"#66c2a5\", colorant\"#fc8d62\", colorant\"#8da0cb\", colorant\"#e78ac3\", \n",
" colorant\"#a6d854\", colorant\"#ffd92f\", colorant\"#e5c494\", colorant\"#b3b3b3\"],\n",
" lchoices=Float64[58, 45, 72.5, 90],\n",
" transform=c -> deuteranopic(c, 0.1),\n",
" cchoices=Float64[20,40],\n",
" hchoices=[75,51,35,120,180,210,270,310]\n",
" )\n",
" convert(Vector{Color}, cs)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import the data from Douven and Schupbach (2015a):"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data = CSV.read(\"data.csv\", DataFrame);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Add difference in explanatory goodness as well as identifiers for participants:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data[!, :diff] = Array{Float64}(data[!, :Judge_A] .- data[!, :Judge_B])\n",
"data[!, :id] = CategoricalArray(repeat(1:26, inner=10))\n",
"first(data, 6)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Below the model that came out on top in Douven and Schupbach's analysis, MMOD:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"form = @formula(Post_AS ~ Post_AO + diff + (Post_AO + diff | id));"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mm = LinearMixedModel(form, data);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fit!(mm)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Below the Bayesian model, MMO:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"form₀ = @formula(Post_AS ~ Post_AO + (Post_AO | id));"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mm₀ = LinearMixedModel(form₀, data);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fit!(mm₀)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, the model MMOAB:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"form₁ = @formula(Post_AS ~ Post_AO + Judge_A + Judge_B + (Post_AO + Judge_A + Judge_B | id));"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mm₁ = LinearMixedModel(form₁, data);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fit!(mm₁)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data1 = copy(data)\n",
"\n",
"ProbEvA = Array{Float64}(undef, 10, 26)\n",
"ProbEvB = Array{Float64}(undef, 10, 26)\n",
"\n",
"for i in 1:26\n",
" \n",
" Pr = vcat(0.5, data1[!, :Post_AO][(i * 10 - 9):(i * 10)])\n",
"\n",
" for j in 1:10\n",
" Pr[j + 1] < Pr[j] ? (ProbEvA[j, i], ProbEvB[j, i]) = (0.25, 0.625) : (ProbEvA[j, i], ProbEvB[j, i]) = (0.75, 0.375)\n",
" end\n",
"\n",
"end\n",
"\n",
"data1[!, :ProbEv] = mean(hcat(cumprod(ProbEvA, dims=1)[:], cumprod(ProbEvB, dims=1)[:]), dims=2)[:]\n",
"data1[!, :HA_Ev] = data1[!, :Post_AO] .* data1[!, :ProbEv];"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"first(data1, 6)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"function aic_vals(d)\n",
" data1[!, :Prob_noise] = ((1-2d)^2 .* data1[!, :HA_Ev] .+ \n",
" (d-2d^2) .* (0.5 .+ data1[!, :ProbEv]) .+ d^2) ./ ((1-2d) .* data1[!, :ProbEv] .+ d)\n",
" BAY = aic(fit!(LinearMixedModel(@formula(Post_AS ~ Prob_noise + (1 + Prob_noise | id)), data1)))\n",
" EXPL = aic(fit!(LinearMixedModel(@formula(Post_AS ~ Prob_noise + Judge_A + Judge_B + \n",
" (1 + Prob_noise + Judge_A + Judge_B | id)), data1)))\n",
" DIFF = aic(fit!(LinearMixedModel(@formula(Post_AS ~ Prob_noise + diff + (1 + Prob_noise + diff | id)), data1)))\n",
" return vcat(BAY, EXPL, DIFF)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"aicVals = Array{Float64, 2}(undef, 100, 3)\n",
"\n",
"for i in 1:100\n",
" aicVals[i,:] = aic_vals(collect(0.0:0.005:0.495)[i])\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"AV = DataFrame(AIC = aicVals[:],\n",
" d_val = repeat(collect(0.000:0.005:0.495), outer=3),\n",
" Predictors = repeat([\"<i>f</i>(O, <i>d</i>)\", \"<i>f</i>(O, <i>d</i>)AB\", \"<i>f</i>(O, <i>d</i>)D\"], inner=100));"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"function bic_vals(d)\n",
" data1[!, :Prob_noise] = ((1-2d)^2 .* data1[!, :HA_Ev] .+ \n",
" (d-2d^2) .* (0.5 .+ data1[!, :ProbEv]) .+ d^2) ./ ((1-2d) .* data1[!, :ProbEv] .+ d)\n",
" BAY = bic(fit!(LinearMixedModel(@formula(Post_AS ~ Prob_noise + (1 + Prob_noise | id)), data1)))\n",
" EXPL = bic(fit!(LinearMixedModel(@formula(Post_AS ~ Prob_noise + Judge_A + Judge_B + \n",
" (1 + Prob_noise + Judge_A + Judge_B | id)), data1)))\n",
" DIFF = bic(fit!(LinearMixedModel(@formula(Post_AS ~ Prob_noise + diff + (1 + Prob_noise + diff | id)), data1)))\n",
" return vcat(BAY, EXPL, DIFF)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"bicVals = Array{Float64, 2}(undef, 100, 3)\n",
"\n",
"for i in 1:100\n",
" bicVals[i,:] = bic_vals(collect(0.0:0.005:0.495)[i])\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"BV = DataFrame(BIC = bicVals[:],\n",
" d_val = repeat(collect(0.000:0.005:0.495), outer=3),\n",
" Predictors = repeat([\"<i>f</i>(O, <i>d</i>)\", \"<i>f</i>(O, <i>d</i>)AB\", \"<i>f</i>(O, <i>d</i>)D\"], inner=100));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Visualize the outcomes:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"CWaic = plot(AV, x=:d_val, y=:AIC, color=:Predictors, Geom.point,\n",
" Guide.ylabel(\"AIC value\"),\n",
" Guide.xlabel(\"<i>d</i>\"),\n",
" Scale.color_discrete_manual(gen_brew_colors(3)...),\n",
" Theme(point_size=2.5pt, minor_label_font_size=10pt,\n",
" major_label_font_size=14pt,\n",
" minor_label_color=colorant\"black\",\n",
" major_label_color=colorant\"black\",\n",
" key_title_font_size=13pt,\n",
" key_label_font_size=11pt,\n",
" key_label_color=colorant\"black\",\n",
" key_title_color=colorant\"black\",\n",
" highlight_width=0.25pt))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"CWbic = plot(BV, x=:d_val, y=:BIC, color=:Predictors, Geom.point,\n",
" Guide.ylabel(\"BIC value\"),\n",
" Guide.xlabel(\"<i>d</i>\"),\n",
" Scale.color_discrete_manual(gen_brew_colors(3)...),\n",
" Theme(point_size=2.5pt, minor_label_font_size=10pt,\n",
" major_label_font_size=14pt,\n",
" minor_label_color=colorant\"black\",\n",
" major_label_color=colorant\"black\",\n",
" key_title_font_size=13pt,\n",
" key_label_font_size=11pt,\n",
" key_label_color=colorant\"black\",\n",
" key_title_color=colorant\"black\",\n",
" highlight_width=0.25pt))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.5.3",
"language": "julia",
"name": "julia-1.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Loading

0 comments on commit 8a6edfb

Please sign in to comment.