diff --git a/.nojekyll b/.nojekyll index 06e546c..d4e09d4 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -027eca6b \ No newline at end of file +0cccd8af \ No newline at end of file diff --git a/search.json b/search.json index f74a939..b532eb6 100644 --- a/search.json +++ b/search.json @@ -1181,7 +1181,7 @@ "href": "sleepstudy.html", "title": "Analysis of the sleepstudy data", "section": "", - "text": "The sleepstudy data are from a study of the effects of sleep deprivation on response time reported in Balkin et al. (2000) and in Belenky et al. (2003). Eighteen subjects were allowed only 3 hours of time to sleep each night for 9 successive nights. Their reaction time was measured each day, starting the day before the first night of sleep deprivation, when the subjects were on their regular sleep schedule.\n\n\n\n\n\n\nNote\n\n\n\nThis description is inaccurate. In fact the first two days were acclimatization, the third was a baseline and sleep deprivation was only enforced after day 2. To allow for comparison with earlier analyses of these data we retain the old data description for this notebook only.\n\n\n\n1 Loading the data\nFirst attach the MixedModels package and other packages for plotting. The CairoMakie package allows the Makie graphics system (Danisch & Krumbiegel, 2021) to generate high quality static images. Activate that package with the SVG (Scalable Vector Graphics) backend.\n\n\nCode\nusing AlgebraOfGraphics\nusing CairoMakie # graphics back-end\nusing DataFrames\nusing KernelDensity # density estimation\nusing MixedModels\nusing MixedModelsMakie # diagnostic plots\nusing Random # random number generators\nusing RCall # call R from Julia\nusing MixedModelsMakie: simplelinreg\nusing SMLP2024: dataset\n\n\nThe sleepstudy data are one of the datasets available with the MixedModels package. It is re-exported by the SMLP2024 package’s dataset function.\n\nsleepstudy = DataFrame(dataset(\"sleepstudy\"))\n\n180×3 DataFrame155 rows omitted\n\n\n\nRow\nsubj\ndays\nreaction\n\n\n\nString\nInt8\nFloat64\n\n\n\n\n1\nS308\n0\n249.56\n\n\n2\nS308\n1\n258.705\n\n\n3\nS308\n2\n250.801\n\n\n4\nS308\n3\n321.44\n\n\n5\nS308\n4\n356.852\n\n\n6\nS308\n5\n414.69\n\n\n7\nS308\n6\n382.204\n\n\n8\nS308\n7\n290.149\n\n\n9\nS308\n8\n430.585\n\n\n10\nS308\n9\n466.353\n\n\n11\nS309\n0\n222.734\n\n\n12\nS309\n1\n205.266\n\n\n13\nS309\n2\n202.978\n\n\n⋮\n⋮\n⋮\n⋮\n\n\n169\nS371\n8\n350.781\n\n\n170\nS371\n9\n369.469\n\n\n171\nS372\n0\n269.412\n\n\n172\nS372\n1\n273.474\n\n\n173\nS372\n2\n297.597\n\n\n174\nS372\n3\n310.632\n\n\n175\nS372\n4\n287.173\n\n\n176\nS372\n5\n329.608\n\n\n177\nS372\n6\n334.482\n\n\n178\nS372\n7\n343.22\n\n\n179\nS372\n8\n369.142\n\n\n180\nS372\n9\n364.124\n\n\n\n\n\n\nFigure 1 displays the data in a multi-panel plot created with the lattice package in R (Sarkar, 2008), using RCall.jl.\n\n\nCode\nlet f = Figure(; size=(800,400))\n yrange = maximum(sleepstudy.reaction) - minimum(sleepstudy.reaction)\n xrange = maximum(sleepstudy.days) - minimum(sleepstudy.days)\n \n reg = combine(groupby(sleepstudy, :subj), \n [:days, :reaction] => NamedTuple{(:intercept, :slope)} ∘ simplelinreg => AsTable)\n sort!(reg, :intercept)\n\n # order of grid positions to plot the facets in\n gridpos = Dict{String, NTuple{2,Int}}()\n for (i, subj) in enumerate(reg.subj)\n gridpos[subj] = fldmod1(i, 9)\n end\n gridpos\n\n axes = Axis[]\n\n # set up all the axes and plot the simple regression lines\n for row in eachrow(reg)\n pos = gridpos[row.subj]\n ax = Axis(f[pos...]; title=row.subj, \n autolimitaspect=xrange/yrange)\n if pos[1] == 1\n hidexdecorations!(ax; grid=false, ticks=false)\n end\n if pos[2] != 1\n hideydecorations!(ax; grid=false, ticks=true)\n end\n push!(axes, ax)\n ablines!(ax, row.intercept, row.slope)\n end\n\n # scatter plot in each facet\n for (grouping, gdf) in pairs(groupby(sleepstudy, :subj))\n pos = gridpos[grouping.subj]\n scatter!(f[pos...], gdf.days, gdf.reaction)\n end\n Label(f[end+1, :], \"Days of sleep deprivation\"; \n tellwidth=false, tellheight=true)\n Label(f[:, 0], \"Average reaction time (ms)\"; \n tellwidth=true, tellheight=false, rotation=pi/2)\n \n linkaxes!(axes...)\n\n # tweak the layout a little\n rowgap!(f.layout, 0)\n colgap!(f.layout, 3)\n colsize!(f.layout, 0, 25)\n rowsize!(f.layout, 1, 100)\n rowsize!(f.layout, 2, 100)\n rowsize!(f.layout, 3, 25)\n f\nend\n\n\n\n\n\n\n\nFigure 1: Average response time versus days of sleep deprivation by subject\n\n\n\n\nEach panel shows the data from one subject and a line fit by least squares to that subject’s data. Starting at the lower left panel and proceeding across rows, the panels are ordered by increasing intercept of the least squares line.\nThere are some deviations from linearity within the panels but the deviations are neither substantial nor systematic.\n\n\n2 Fitting an initial model\n\ncontrasts = Dict{Symbol,Any}(:subj => Grouping())\nm1 = let f = @formula(reaction ~ 1 + days + (1 + days | subj))\n fit(MixedModel, f, sleepstudy; contrasts)\nend\n\n\n\n\n\nEst.\nSE\nz\np\nσ_subj\n\n\n(Intercept)\n251.4051\n6.6323\n37.91\n<1e-99\n23.7805\n\n\ndays\n10.4673\n1.5022\n6.97\n<1e-11\n5.7168\n\n\nResidual\n25.5918\n\n\n\n\n\n\n\n\n\nThis model includes fixed effects for the intercept, representing the typical reaction time at the beginning of the experiment with zero days of sleep deprivation, and the slope w.r.t. days of sleep deprivation. The parameter estimates are about 250 ms. typical reaction time without deprivation and a typical increase of 10.5 ms. per day of sleep deprivation.\nThe random effects represent shifts from the typical behavior for each subject. The shift in the intercept has a standard deviation of about 24 ms. which would suggest a range of about 200 ms. to 300 ms. in the intercepts. Similarly within-subject slopes would be expected to have a range of about 0 ms./day up to 20 ms./day.\nThe random effects for the slope and for the intercept are allowed to be correlated within subject. The estimated correlation, 0.08, is small. This estimate is not shown in the default display above but is shown in the output from VarCorr (variance components and correlations).\n\nVarCorr(m1)\n\n\n\n\n\nColumn\nVariance\nStd.Dev\nCorr.\n\n\nsubj\n(Intercept)\n565.51065\n23.78047\n\n\n\n\ndays\n32.68212\n5.71683\n+0.08\n\n\nResidual\n\n654.94145\n25.59182\n\n\n\n\n\n\nTechnically, the random effects for each subject are unobserved random variables and are not “parameters” in the model per se. Hence we do not report standard errors or confidence intervals for these deviations. However, we can produce prediction intervals on the random effects for each subject. Because the experimental design is balanced, these intervals will have the same width for all subjects.\nA plot of the prediction intervals versus the level of the grouping factor (subj, in this case) is sometimes called a caterpillar plot because it can look like a fuzzy caterpillar if there are many levels of the grouping factor. By default, the levels of the grouping factor are sorted by increasing value of the first random effect.\n\n\nCode\ncaterpillar(m1; vline_at_zero=true)\n\n\n\n\n\n\n\nFigure 2: Prediction intervals on random effects for model m1\n\n\n\n\nFigure 2 reinforces the conclusion that there is little correlation between the random effect for intercept and the random effect for slope.\n\n\n3 A model with uncorrelated random effects\nThe zerocorr function applied to a random-effects term creates uncorrelated vector-valued per-subject random effects.\n\nm2 = let f = @formula reaction ~ 1 + days + zerocorr(1 + days | subj)\n fit(MixedModel, f, sleepstudy; contrasts)\nend\n\n\n\n\n\nEst.\nSE\nz\np\nσ_subj\n\n\n(Intercept)\n251.4051\n6.7077\n37.48\n<1e-99\n24.1714\n\n\ndays\n10.4673\n1.5193\n6.89\n<1e-11\n5.7994\n\n\nResidual\n25.5561\n\n\n\n\n\n\n\n\n\nAgain, the default display doesn’t show that there is no correlation parameter to be estimated in this model, but the VarCorr display does.\n\nVarCorr(m2)\n\n\n\n\n\nColumn\nVariance\nStd.Dev\nCorr.\n\n\nsubj\n(Intercept)\n584.25897\n24.17145\n\n\n\n\ndays\n33.63281\n5.79938\n.\n\n\nResidual\n\n653.11578\n25.55613\n\n\n\n\n\n\nThis model has a slightly lower log-likelihood than does m1 and one fewer parameter than m1. A likelihood-ratio test can be used to compare these nested models.\n\nMixedModels.likelihoodratiotest(m2, m1)\n\n\n\n\n\nmodel-dof\ndeviance\nχ²\nχ²-dof\nP(>χ²)\n\n\nreaction ~ 1 + days + zerocorr(1 + days | subj)\n5\n1752\n\n\n\n\n\nreaction ~ 1 + days + (1 + days | subj)\n6\n1752\n0\n1\n0.8004\n\n\n\n\n\nAlternatively, the AIC or BIC values can be compared.\n\n\nCode\nlet mods = [m2, m1]\n Table(;\n model=[:m2, :m1],\n pars=dof.(mods),\n geomdof=(sum ∘ leverage).(mods),\n AIC=aic.(mods),\n BIC=bic.(mods),\n AICc=aicc.(mods),\n )\nend\n\n\nTable with 6 columns and 2 rows:\n model pars geomdof AIC BIC AICc\n ┌────────────────────────────────────────────────\n 1 │ m2 5 29.045 1762.0 1777.97 1762.35\n 2 │ m1 6 28.6115 1763.94 1783.1 1764.42\n\n\nThe goodness of fit measures: AIC, BIC, and AICc, are all on a “smaller is better” scale and, hence, they all prefer m2.\nThe pars column, which is the same as the model-dof column in the likelihood ratio test output, is simply a count of the number of parameters to be estimated when fitting the model. For example, in m2 there are two fixed-effects parameters and three variance components (including the residual variance).\nAn alternative, more geometrically inspired definition of “degrees of freedom”, is the sum of the leverage values, called geomdof in this table.\nInterestingly, the model with fewer parameters, m2, has a greater sum of the leverage values than the model with more parameters, m1. We’re not sure what to make of that.\nIn both cases the sum of the leverage values is toward the upper end of the range of possible values, which is the rank of the fixed-effects model matrix (2) up to the rank of the fixed-effects plus the random effects model matrix (2 + 36 = 38).\n\n\n\n\n\n\nNote\n\n\n\nI think that the upper bound may be 36, not 38, because the two columns of X lie in the column span of Z\n\n\nThis comparison does show, however, that a simple count of the parameters in a mixed-effects model can underestimate, sometimes drastically, the model complexity. This is because a single variance component or multiple components can add many dimensions to the linear predictor.\n\n\n4 Some diagnostic plots\nIn mixed-effects models the linear predictor expression incorporates fixed-effects parameters, which summarize trends for the population or certain well-defined subpopulations, and random effects which represent deviations associated with the experimental units or observational units - individual subjects, in this case. The random effects are modeled as unobserved random variables.\nThe conditional means of these random variables, sometimes called the BLUPs or Best Linear Unbiased Predictors, are not simply the least squares estimates. They are attenuated or shrunk towards zero to reflect the fact that the individuals are assumed to come from a population. A shrinkage plot, Figure 3, shows the BLUPs from the model fit compared to the values without any shrinkage. If the BLUPs are similar to the unshrunk values then the more complicated model accounting for individual differences is supported. If the BLUPs are strongly shrunk towards zero then the additional complexity in the model to account for individual differences is not providing sufficient increase in fidelity to the data to warrant inclusion.\n\n\nCode\nshrinkageplot!(Figure(; resolution=(500, 500)), m1)\n\n\n┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.\n└ @ Makie ~/.julia/packages/Makie/WgbrE/src/scenes.jl:227\n\n\n\n\n\n\n\nFigure 3: Shrinkage plot of means of the random effects in model m1\n\n\n\n\n\n\n\n\n\n\nNote\n\n\n\nThis plot could be drawn as shrinkageplot(m1). The reason for explicitly creating a Figure to be modified by shrinkageplot! is to control the resolution.\n\n\nThis plot shows an intermediate pattern. The random effects are somewhat shrunk toward the origin, a model simplification trend, but not completely shrunk - indicating that fidelity to the data is enhanced with these additional coefficients in the linear predictor.\nIf the shrinkage were primarily in one direction - for example, if the arrows from the unshrunk values to the shrunk values were mostly in the vertical direction - then we would get an indication that we could drop the random effect for slope and revert to a simpler model. This is not the case here.\nAs would be expected, the unshrunk values that are further from the origin tend to be shrunk more toward the origin. That is, the arrows that originate furthest from the origin are longer. However, that is not always the case. The arrow in the upper right corner, from S337, is relatively short. Examination of the panel for S337 in the data plot shows a strong linear trend, even though both the intercept and the slope are unusually large. The neighboring panels in the data plot, S330 and S331, have more variability around the least squares line and are subject to a greater amount of shrinkage in the model. (They correspond to the two arrows on the right hand side of the figure around -5 on the vertical scale.)\n\n\n5 Assessing variability by bootstrapping\nThe speed of fitting linear mixed-effects models using MixedModels.jl allows for using simulation-based approaches to inference instead of relying on approximate standard errors. A parametric bootstrap sample for model m is a collection of models of the same form as m fit to data values simulated from m. That is, we pretend that m and its parameter values are the true parameter values, simulate data from these values, and estimate parameters from the simulated data.\nSimulating and fitting a substantial number of model fits, 5000 in this case, takes only a few seconds, following which we extract a data frame of the parameter estimates and plot densities of some of these estimates.\n\nrng = Random.seed!(42) # initialize a random number generator\nm1bstp = parametricbootstrap(rng, 5000, m1)\ntbl = m1bstp.tbl\n\nTable with 10 columns and 5000 rows:\n obj β1 β2 σ σ1 σ2 ρ1 ⋯\n ┌────────────────────────────────────────────────────────────────────\n 1 │ 1717.29 260.712 9.84975 23.4092 15.3314 6.40292 -0.0259483 ⋯\n 2 │ 1744.06 262.253 12.3008 25.7047 16.3186 5.54691 0.552569 ⋯\n 3 │ 1714.16 253.149 12.879 22.2753 25.4787 6.1444 0.0691544 ⋯\n 4 │ 1711.54 263.376 11.5798 23.3128 18.8039 4.65569 0.103361 ⋯\n 5 │ 1741.66 248.429 9.39444 25.4355 20.1412 5.27358 -0.163606 ⋯\n 6 │ 1754.81 256.794 8.024 26.5088 10.6781 7.14155 0.335447 ⋯\n 7 │ 1777.73 253.388 8.83556 27.8623 17.8325 7.17386 0.00379005 ⋯\n 8 │ 1768.59 254.441 11.4479 27.4033 16.2484 6.67057 0.725411 ⋯\n 9 │ 1753.56 244.906 11.3423 25.6046 25.3607 5.98654 -0.171821 ⋯\n 10 │ 1722.61 257.088 9.18397 23.3386 24.9274 5.18012 0.181143 ⋯\n 11 │ 1738.16 251.262 11.6568 25.7823 17.6663 4.07215 0.258142 ⋯\n 12 │ 1747.76 258.302 12.8015 26.1083 19.2395 5.06103 0.879669 ⋯\n 13 │ 1745.91 254.57 11.8062 24.8863 24.2513 6.14642 0.0126995 ⋯\n 14 │ 1738.8 251.179 10.3226 24.2672 23.7195 6.32645 0.368592 ⋯\n 15 │ 1724.76 238.603 11.5045 25.2301 19.026 3.64035 -0.346558 ⋯\n 16 │ 1777.7 254.133 8.26397 26.9846 26.3715 7.8283 -0.288773 ⋯\n 17 │ 1748.33 251.571 9.5294 26.2927 21.9611 4.31316 -0.150104 ⋯\n ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱\n\n\nAn empirical density plot of the estimates for the fixed-effects coefficients, Figure 4, shows the normal distribution, “bell-curve”, shape as we might expect.\n\n\nCode\nbegin\n f1 = Figure(; resolution=(1000, 400))\n CairoMakie.density!(\n Axis(f1[1, 1]; xlabel=\"Intercept [ms]\"), tbl.β1\n )\n CairoMakie.density!(\n Axis(f1[1, 2]; xlabel=\"Coefficient of days [ms/day]\"),\n tbl.β2\n )\n f1\nend\n\n\n┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.\n└ @ Makie ~/.julia/packages/Makie/WgbrE/src/scenes.jl:227\n\n\n\n\n\n\n\nFigure 4: Empirical density plots of bootstrap replications of fixed-effects parameter estimates\n\n\n\n\nIt is also possible to create interval estimates of the parameters from the bootstrap replicates. We define the 1-α shortestcovint to be the shortest interval that contains a proportion 1-α (defaults to 95%) of the bootstrap estimates of the parameter.\n\nTable(shortestcovint(m1bstp))\n\nTable with 5 columns and 6 rows:\n type group names lower upper\n ┌─────────────────────────────────────────────────────\n 1 │ β missing (Intercept) 239.64 265.228\n 2 │ β missing days 7.42347 13.1607\n 3 │ σ subj (Intercept) 10.1722 33.0876\n 4 │ σ subj days 2.99474 7.66118\n 5 │ ρ subj (Intercept), days -0.40135 1.0\n 6 │ σ residual missing 22.701 28.5016\n\n\nThe intervals look reasonable except that the upper end point of the interval for ρ1, the correlation coefficient, is 1.0 . It turns out that the estimates of ρ have a great deal of variability.\nBecause there are several values on the boundary (ρ = 1.0) and a pulse like this is not handled well by a density plot, we plot this sample as a histogram, Figure 5.\n\n\nCode\nhist(\n tbl.ρ1;\n bins=40,\n axis=(; xlabel=\"Estimated correlation of the random effects\"),\n figure=(; resolution=(500, 500)),\n)\n\n\n┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.\n└ @ Makie ~/.julia/packages/Makie/WgbrE/src/scenes.jl:227\n\n\n\n\n\n\n\nFigure 5: Histogram of bootstrap replications of the within-subject correlation parameter\n\n\n\n\nFinally, density plots for the variance components (but on the scale of the standard deviation), Figure 6, show reasonable symmetry.\n\n\nCode\nbegin\n f2 = Figure(; resolution=(1000, 300))\n CairoMakie.density!(\n Axis(f2[1, 1]; xlabel=\"Residual σ\"),\n tbl.σ,\n )\n CairoMakie.density!(\n Axis(f2[1, 2]; xlabel=\"subj-Intercept σ\"),\n tbl.σ1,\n )\n CairoMakie.density!(\n Axis(f2[1, 3]; xlabel=\"subj-slope σ\"),\n tbl.σ2,\n )\n f2\nend\n\n\n┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.\n└ @ Makie ~/.julia/packages/Makie/WgbrE/src/scenes.jl:227\n\n\n\n\n\n\n\nFigure 6: Empirical density plots of bootstrap replicates of standard deviation estimates\n\n\n\n\nThe estimates of the coefficients, β₁ and β₂, are not highly correlated as shown in a scatterplot of the bootstrap estimates, Figure 7 .\n\nvcov(m1; corr=true) # correlation estimate from the model\n\n2×2 Matrix{Float64}:\n 1.0 -0.137545\n -0.137545 1.0\n\n\n\n\nCode\nlet\n scatter(\n tbl.β1, tbl.β2,\n color=(:blue, 0.20),\n axis=(; xlabel=\"Intercept\", ylabel=\"Coefficient of days\"),\n figure=(; resolution=(500, 500)),\n )\n contour!(kde((tbl.β1, tbl.β2)))\n current_figure()\nend\n\n\n┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.\n└ @ Makie ~/.julia/packages/Makie/WgbrE/src/scenes.jl:227\n\n\n\n\n\n\n\nFigure 7: Scatter-plot of bootstrap replicates of fixed-effects estimates with contours\n\n\n\n\n\n\n6 References\n\n\nBalkin, T., Thome, D., Sing, H., Thomas, M., Redmond, D., Wesensten, N., Williams, J., Hall, S., & Belenky, G. (2000). Effects of sleep schedules on commercial motor vehicle driver performance (DOT-MC-00-133). Federal Motor Carrier Safety Administration. https://doi.org/10.21949/1503015.\n\n\nBelenky, G., Wesensten, N. J., Thorne, D. R., Thomas, M. L., Sing, H. C., Redmond, D. P., Russo, M. B., & Balkin, T. J. (2003). Patterns of performance degradation and restoration during sleep restriction and subsequent recovery: A sleep dose-response study. Journal of Sleep Research, 12(1), 1–12. https://doi.org/10.1046/j.1365-2869.2003.00337.x\n\n\nDanisch, S., & Krumbiegel, J. (2021). Makie.jl: Flexible high-performance data visualization for julia. Journal of Open Source Software, 6(65), 3349. https://doi.org/10.21105/joss.03349\n\n\nSarkar, D. (2008). Lattice: Mutivariate data visualization with r. Springer-Verlag GmbH. https://www.ebook.de/de/product/11429038/deepayan_sarkar_lattice.html\n\n\n\n\n\n\n Back to top", + "text": "The sleepstudy data are from a study of the effects of sleep deprivation on response time reported in Balkin et al. (2000) and in Belenky et al. (2003). Eighteen subjects were allowed only 3 hours of time to sleep each night for 9 successive nights. Their reaction time was measured each day, starting the day before the first night of sleep deprivation, when the subjects were on their regular sleep schedule.\n\n\n\n\n\n\nNote\n\n\n\nThis description is inaccurate. In fact the first two days were acclimatization, the third was a baseline and sleep deprivation was only enforced after day 2. To allow for comparison with earlier analyses of these data we retain the old data description for this notebook only.\n\n\n\n1 Loading the data\nFirst attach the MixedModels package and other packages for plotting. The CairoMakie package allows the Makie graphics system (Danisch & Krumbiegel, 2021) to generate high quality static images. Activate that package with the SVG (Scalable Vector Graphics) backend.\n\n\nCode\nusing AlgebraOfGraphics\nusing CairoMakie # graphics back-end\nusing DataFrames\nusing KernelDensity # density estimation\nusing MixedModels\nusing MixedModelsMakie # diagnostic plots\nusing Random # random number generators\nusing RCall # call R from Julia\nusing MixedModelsMakie: simplelinreg\nusing SMLP2024: dataset\n\n\nThe sleepstudy data are one of the datasets available with the MixedModels package. It is re-exported by the SMLP2024 package’s dataset function.\n\nsleepstudy = DataFrame(dataset(\"sleepstudy\"))\n\n180×3 DataFrame155 rows omitted\n\n\n\nRow\nsubj\ndays\nreaction\n\n\n\nString\nInt8\nFloat64\n\n\n\n\n1\nS308\n0\n249.56\n\n\n2\nS308\n1\n258.705\n\n\n3\nS308\n2\n250.801\n\n\n4\nS308\n3\n321.44\n\n\n5\nS308\n4\n356.852\n\n\n6\nS308\n5\n414.69\n\n\n7\nS308\n6\n382.204\n\n\n8\nS308\n7\n290.149\n\n\n9\nS308\n8\n430.585\n\n\n10\nS308\n9\n466.353\n\n\n11\nS309\n0\n222.734\n\n\n12\nS309\n1\n205.266\n\n\n13\nS309\n2\n202.978\n\n\n⋮\n⋮\n⋮\n⋮\n\n\n169\nS371\n8\n350.781\n\n\n170\nS371\n9\n369.469\n\n\n171\nS372\n0\n269.412\n\n\n172\nS372\n1\n273.474\n\n\n173\nS372\n2\n297.597\n\n\n174\nS372\n3\n310.632\n\n\n175\nS372\n4\n287.173\n\n\n176\nS372\n5\n329.608\n\n\n177\nS372\n6\n334.482\n\n\n178\nS372\n7\n343.22\n\n\n179\nS372\n8\n369.142\n\n\n180\nS372\n9\n364.124\n\n\n\n\n\n\nFigure 1 displays the data in a multi-panel plot created with the lattice package in R (Sarkar, 2008), using RCall.jl.\n\n\nCode\nlet f = Figure(; size=(700, 400))\n yrange = maximum(sleepstudy.reaction) - minimum(sleepstudy.reaction)\n xrange = maximum(sleepstudy.days) - minimum(sleepstudy.days)\n \n reg = combine(groupby(sleepstudy, :subj), \n [:days, :reaction] => NamedTuple{(:intercept, :slope)} ∘ simplelinreg => AsTable)\n sort!(reg, :intercept)\n\n # order of grid positions to plot the facets in\n gridpos = Dict{String, NTuple{2,Int}}()\n for (i, subj) in enumerate(reg.subj)\n gridpos[subj] = fldmod1(i, 9)\n end\n gridpos\n\n axes = Axis[]\n\n # set up all the axes and plot the simple regression lines\n for row in eachrow(reg)\n pos = gridpos[row.subj]\n ax = Axis(f[pos...]; title=row.subj, \n autolimitaspect=xrange/yrange)\n if pos[1] == 1\n hidexdecorations!(ax; grid=false, ticks=false)\n end\n if pos[2] != 1\n hideydecorations!(ax; grid=false, ticks=true)\n end\n push!(axes, ax)\n ablines!(ax, row.intercept, row.slope)\n end\n\n # scatter plot in each facet\n for (grouping, gdf) in pairs(groupby(sleepstudy, :subj))\n pos = gridpos[grouping.subj]\n scatter!(f[pos...], gdf.days, gdf.reaction)\n end\n Label(f[end+1, :], \"Days of sleep deprivation\"; \n tellwidth=false, tellheight=true)\n Label(f[:, 0], \"Average reaction time (ms)\"; \n tellwidth=true, tellheight=false, rotation=pi/2)\n \n linkaxes!(axes...)\n\n # tweak the layout a little\n rowgap!(f.layout, 0)\n colgap!(f.layout, 3)\n colsize!(f.layout, 0, 25)\n rowsize!(f.layout, 1, 100)\n rowsize!(f.layout, 2, 100)\n rowsize!(f.layout, 3, 25)\n f\nend\n\n\n\n\n\n\n\nFigure 1: Average response time versus days of sleep deprivation by subject\n\n\n\n\nEach panel shows the data from one subject and a line fit by least squares to that subject’s data. Starting at the lower left panel and proceeding across rows, the panels are ordered by increasing intercept of the least squares line.\nThere are some deviations from linearity within the panels but the deviations are neither substantial nor systematic.\n\n\n2 Fitting an initial model\n\ncontrasts = Dict{Symbol,Any}(:subj => Grouping())\nm1 = let f = @formula(reaction ~ 1 + days + (1 + days | subj))\n fit(MixedModel, f, sleepstudy; contrasts)\nend\n\n\n\n\n\nEst.\nSE\nz\np\nσ_subj\n\n\n(Intercept)\n251.4051\n6.6323\n37.91\n<1e-99\n23.7805\n\n\ndays\n10.4673\n1.5022\n6.97\n<1e-11\n5.7168\n\n\nResidual\n25.5918\n\n\n\n\n\n\n\n\n\nThis model includes fixed effects for the intercept, representing the typical reaction time at the beginning of the experiment with zero days of sleep deprivation, and the slope w.r.t. days of sleep deprivation. The parameter estimates are about 250 ms. typical reaction time without deprivation and a typical increase of 10.5 ms. per day of sleep deprivation.\nThe random effects represent shifts from the typical behavior for each subject. The shift in the intercept has a standard deviation of about 24 ms. which would suggest a range of about 200 ms. to 300 ms. in the intercepts. Similarly within-subject slopes would be expected to have a range of about 0 ms./day up to 20 ms./day.\nThe random effects for the slope and for the intercept are allowed to be correlated within subject. The estimated correlation, 0.08, is small. This estimate is not shown in the default display above but is shown in the output from VarCorr (variance components and correlations).\n\nVarCorr(m1)\n\n\n\n\n\nColumn\nVariance\nStd.Dev\nCorr.\n\n\nsubj\n(Intercept)\n565.51065\n23.78047\n\n\n\n\ndays\n32.68212\n5.71683\n+0.08\n\n\nResidual\n\n654.94145\n25.59182\n\n\n\n\n\n\nTechnically, the random effects for each subject are unobserved random variables and are not “parameters” in the model per se. Hence we do not report standard errors or confidence intervals for these deviations. However, we can produce prediction intervals on the random effects for each subject. Because the experimental design is balanced, these intervals will have the same width for all subjects.\nA plot of the prediction intervals versus the level of the grouping factor (subj, in this case) is sometimes called a caterpillar plot because it can look like a fuzzy caterpillar if there are many levels of the grouping factor. By default, the levels of the grouping factor are sorted by increasing value of the first random effect.\n\n\nCode\ncaterpillar(m1; vline_at_zero=true)\n\n\n\n\n\n\n\nFigure 2: Prediction intervals on random effects for model m1\n\n\n\n\nFigure 2 reinforces the conclusion that there is little correlation between the random effect for intercept and the random effect for slope.\n\n\n3 A model with uncorrelated random effects\nThe zerocorr function applied to a random-effects term creates uncorrelated vector-valued per-subject random effects.\n\nm2 = let f = @formula reaction ~ 1 + days + zerocorr(1 + days | subj)\n fit(MixedModel, f, sleepstudy; contrasts)\nend\n\n\n\n\n\nEst.\nSE\nz\np\nσ_subj\n\n\n(Intercept)\n251.4051\n6.7077\n37.48\n<1e-99\n24.1714\n\n\ndays\n10.4673\n1.5193\n6.89\n<1e-11\n5.7994\n\n\nResidual\n25.5561\n\n\n\n\n\n\n\n\n\nAgain, the default display doesn’t show that there is no correlation parameter to be estimated in this model, but the VarCorr display does.\n\nVarCorr(m2)\n\n\n\n\n\nColumn\nVariance\nStd.Dev\nCorr.\n\n\nsubj\n(Intercept)\n584.25897\n24.17145\n\n\n\n\ndays\n33.63281\n5.79938\n.\n\n\nResidual\n\n653.11578\n25.55613\n\n\n\n\n\n\nThis model has a slightly lower log-likelihood than does m1 and one fewer parameter than m1. A likelihood-ratio test can be used to compare these nested models.\n\nMixedModels.likelihoodratiotest(m2, m1)\n\n\n\n\n\nmodel-dof\ndeviance\nχ²\nχ²-dof\nP(>χ²)\n\n\nreaction ~ 1 + days + zerocorr(1 + days | subj)\n5\n1752\n\n\n\n\n\nreaction ~ 1 + days + (1 + days | subj)\n6\n1752\n0\n1\n0.8004\n\n\n\n\n\nAlternatively, the AIC or BIC values can be compared.\n\n\nCode\nlet mods = [m2, m1]\n Table(;\n model=[:m2, :m1],\n pars=dof.(mods),\n geomdof=(sum ∘ leverage).(mods),\n AIC=aic.(mods),\n BIC=bic.(mods),\n AICc=aicc.(mods),\n )\nend\n\n\nTable with 6 columns and 2 rows:\n model pars geomdof AIC BIC AICc\n ┌────────────────────────────────────────────────\n 1 │ m2 5 29.045 1762.0 1777.97 1762.35\n 2 │ m1 6 28.6115 1763.94 1783.1 1764.42\n\n\nThe goodness of fit measures: AIC, BIC, and AICc, are all on a “smaller is better” scale and, hence, they all prefer m2.\nThe pars column, which is the same as the model-dof column in the likelihood ratio test output, is simply a count of the number of parameters to be estimated when fitting the model. For example, in m2 there are two fixed-effects parameters and three variance components (including the residual variance).\nAn alternative, more geometrically inspired definition of “degrees of freedom”, is the sum of the leverage values, called geomdof in this table.\nInterestingly, the model with fewer parameters, m2, has a greater sum of the leverage values than the model with more parameters, m1. We’re not sure what to make of that.\nIn both cases the sum of the leverage values is toward the upper end of the range of possible values, which is the rank of the fixed-effects model matrix (2) up to the rank of the fixed-effects plus the random effects model matrix (2 + 36 = 38).\n\n\n\n\n\n\nNote\n\n\n\nI think that the upper bound may be 36, not 38, because the two columns of X lie in the column span of Z\n\n\nThis comparison does show, however, that a simple count of the parameters in a mixed-effects model can underestimate, sometimes drastically, the model complexity. This is because a single variance component or multiple components can add many dimensions to the linear predictor.\n\n\n4 Some diagnostic plots\nIn mixed-effects models the linear predictor expression incorporates fixed-effects parameters, which summarize trends for the population or certain well-defined subpopulations, and random effects which represent deviations associated with the experimental units or observational units - individual subjects, in this case. The random effects are modeled as unobserved random variables.\nThe conditional means of these random variables, sometimes called the BLUPs or Best Linear Unbiased Predictors, are not simply the least squares estimates. They are attenuated or shrunk towards zero to reflect the fact that the individuals are assumed to come from a population. A shrinkage plot, Figure 3, shows the BLUPs from the model fit compared to the values without any shrinkage. If the BLUPs are similar to the unshrunk values then the more complicated model accounting for individual differences is supported. If the BLUPs are strongly shrunk towards zero then the additional complexity in the model to account for individual differences is not providing sufficient increase in fidelity to the data to warrant inclusion.\n\n\nCode\nshrinkageplot!(Figure(; resolution=(500, 500)), m1)\n\n\n┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.\n└ @ Makie ~/.julia/packages/Makie/WgbrE/src/scenes.jl:227\n\n\n\n\n\n\n\nFigure 3: Shrinkage plot of means of the random effects in model m1\n\n\n\n\n\n\n\n\n\n\nNote\n\n\n\nThis plot could be drawn as shrinkageplot(m1). The reason for explicitly creating a Figure to be modified by shrinkageplot! is to control the resolution.\n\n\nThis plot shows an intermediate pattern. The random effects are somewhat shrunk toward the origin, a model simplification trend, but not completely shrunk - indicating that fidelity to the data is enhanced with these additional coefficients in the linear predictor.\nIf the shrinkage were primarily in one direction - for example, if the arrows from the unshrunk values to the shrunk values were mostly in the vertical direction - then we would get an indication that we could drop the random effect for slope and revert to a simpler model. This is not the case here.\nAs would be expected, the unshrunk values that are further from the origin tend to be shrunk more toward the origin. That is, the arrows that originate furthest from the origin are longer. However, that is not always the case. The arrow in the upper right corner, from S337, is relatively short. Examination of the panel for S337 in the data plot shows a strong linear trend, even though both the intercept and the slope are unusually large. The neighboring panels in the data plot, S330 and S331, have more variability around the least squares line and are subject to a greater amount of shrinkage in the model. (They correspond to the two arrows on the right hand side of the figure around -5 on the vertical scale.)\n\n\n5 Assessing variability by bootstrapping\nThe speed of fitting linear mixed-effects models using MixedModels.jl allows for using simulation-based approaches to inference instead of relying on approximate standard errors. A parametric bootstrap sample for model m is a collection of models of the same form as m fit to data values simulated from m. That is, we pretend that m and its parameter values are the true parameter values, simulate data from these values, and estimate parameters from the simulated data.\nSimulating and fitting a substantial number of model fits, 5000 in this case, takes only a few seconds, following which we extract a data frame of the parameter estimates and plot densities of some of these estimates.\n\nrng = Random.seed!(42) # initialize a random number generator\nm1bstp = parametricbootstrap(rng, 5000, m1)\ntbl = m1bstp.tbl\n\nTable with 10 columns and 5000 rows:\n obj β1 β2 σ σ1 σ2 ρ1 ⋯\n ┌────────────────────────────────────────────────────────────────────\n 1 │ 1717.29 260.712 9.84975 23.4092 15.3314 6.40292 -0.0259483 ⋯\n 2 │ 1744.06 262.253 12.3008 25.7047 16.3186 5.54691 0.552569 ⋯\n 3 │ 1714.16 253.149 12.879 22.2753 25.4787 6.1444 0.0691544 ⋯\n 4 │ 1711.54 263.376 11.5798 23.3128 18.8039 4.65569 0.103361 ⋯\n 5 │ 1741.66 248.429 9.39444 25.4355 20.1412 5.27358 -0.163606 ⋯\n 6 │ 1754.81 256.794 8.024 26.5088 10.6781 7.14155 0.335447 ⋯\n 7 │ 1777.73 253.388 8.83556 27.8623 17.8325 7.17386 0.00379005 ⋯\n 8 │ 1768.59 254.441 11.4479 27.4033 16.2484 6.67057 0.725411 ⋯\n 9 │ 1753.56 244.906 11.3423 25.6046 25.3607 5.98654 -0.171821 ⋯\n 10 │ 1722.61 257.088 9.18397 23.3386 24.9274 5.18012 0.181143 ⋯\n 11 │ 1738.16 251.262 11.6568 25.7823 17.6663 4.07215 0.258142 ⋯\n 12 │ 1747.76 258.302 12.8015 26.1083 19.2395 5.06103 0.879669 ⋯\n 13 │ 1745.91 254.57 11.8062 24.8863 24.2513 6.14642 0.0126995 ⋯\n 14 │ 1738.8 251.179 10.3226 24.2672 23.7195 6.32645 0.368592 ⋯\n 15 │ 1724.76 238.603 11.5045 25.2301 19.026 3.64035 -0.346558 ⋯\n 16 │ 1777.7 254.133 8.26397 26.9846 26.3715 7.8283 -0.288773 ⋯\n 17 │ 1748.33 251.571 9.5294 26.2927 21.9611 4.31316 -0.150104 ⋯\n ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱\n\n\nAn empirical density plot of the estimates for the fixed-effects coefficients, Figure 4, shows the normal distribution, “bell-curve”, shape as we might expect.\n\n\nCode\nbegin\n f1 = Figure(; resolution=(1000, 400))\n CairoMakie.density!(\n Axis(f1[1, 1]; xlabel=\"Intercept [ms]\"), tbl.β1\n )\n CairoMakie.density!(\n Axis(f1[1, 2]; xlabel=\"Coefficient of days [ms/day]\"),\n tbl.β2\n )\n f1\nend\n\n\n┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.\n└ @ Makie ~/.julia/packages/Makie/WgbrE/src/scenes.jl:227\n\n\n\n\n\n\n\nFigure 4: Empirical density plots of bootstrap replications of fixed-effects parameter estimates\n\n\n\n\nIt is also possible to create interval estimates of the parameters from the bootstrap replicates. We define the 1-α shortestcovint to be the shortest interval that contains a proportion 1-α (defaults to 95%) of the bootstrap estimates of the parameter.\n\nTable(shortestcovint(m1bstp))\n\nTable with 5 columns and 6 rows:\n type group names lower upper\n ┌─────────────────────────────────────────────────────\n 1 │ β missing (Intercept) 239.64 265.228\n 2 │ β missing days 7.42347 13.1607\n 3 │ σ subj (Intercept) 10.1722 33.0876\n 4 │ σ subj days 2.99474 7.66118\n 5 │ ρ subj (Intercept), days -0.40135 1.0\n 6 │ σ residual missing 22.701 28.5016\n\n\nThe intervals look reasonable except that the upper end point of the interval for ρ1, the correlation coefficient, is 1.0 . It turns out that the estimates of ρ have a great deal of variability.\nBecause there are several values on the boundary (ρ = 1.0) and a pulse like this is not handled well by a density plot, we plot this sample as a histogram, Figure 5.\n\n\nCode\nhist(\n tbl.ρ1;\n bins=40,\n axis=(; xlabel=\"Estimated correlation of the random effects\"),\n figure=(; resolution=(500, 500)),\n)\n\n\n┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.\n└ @ Makie ~/.julia/packages/Makie/WgbrE/src/scenes.jl:227\n\n\n\n\n\n\n\nFigure 5: Histogram of bootstrap replications of the within-subject correlation parameter\n\n\n\n\nFinally, density plots for the variance components (but on the scale of the standard deviation), Figure 6, show reasonable symmetry.\n\n\nCode\nbegin\n f2 = Figure(; resolution=(1000, 300))\n CairoMakie.density!(\n Axis(f2[1, 1]; xlabel=\"Residual σ\"),\n tbl.σ,\n )\n CairoMakie.density!(\n Axis(f2[1, 2]; xlabel=\"subj-Intercept σ\"),\n tbl.σ1,\n )\n CairoMakie.density!(\n Axis(f2[1, 3]; xlabel=\"subj-slope σ\"),\n tbl.σ2,\n )\n f2\nend\n\n\n┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.\n└ @ Makie ~/.julia/packages/Makie/WgbrE/src/scenes.jl:227\n\n\n\n\n\n\n\nFigure 6: Empirical density plots of bootstrap replicates of standard deviation estimates\n\n\n\n\nThe estimates of the coefficients, β₁ and β₂, are not highly correlated as shown in a scatterplot of the bootstrap estimates, Figure 7 .\n\nvcov(m1; corr=true) # correlation estimate from the model\n\n2×2 Matrix{Float64}:\n 1.0 -0.137545\n -0.137545 1.0\n\n\n\n\nCode\nlet\n scatter(\n tbl.β1, tbl.β2,\n color=(:blue, 0.20),\n axis=(; xlabel=\"Intercept\", ylabel=\"Coefficient of days\"),\n figure=(; resolution=(500, 500)),\n )\n contour!(kde((tbl.β1, tbl.β2)))\n current_figure()\nend\n\n\n┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.\n└ @ Makie ~/.julia/packages/Makie/WgbrE/src/scenes.jl:227\n\n\n\n\n\n\n\nFigure 7: Scatter-plot of bootstrap replicates of fixed-effects estimates with contours\n\n\n\n\n\n\n6 References\n\n\nBalkin, T., Thome, D., Sing, H., Thomas, M., Redmond, D., Wesensten, N., Williams, J., Hall, S., & Belenky, G. (2000). Effects of sleep schedules on commercial motor vehicle driver performance (DOT-MC-00-133). Federal Motor Carrier Safety Administration. https://doi.org/10.21949/1503015.\n\n\nBelenky, G., Wesensten, N. J., Thorne, D. R., Thomas, M. L., Sing, H. C., Redmond, D. P., Russo, M. B., & Balkin, T. J. (2003). Patterns of performance degradation and restoration during sleep restriction and subsequent recovery: A sleep dose-response study. Journal of Sleep Research, 12(1), 1–12. https://doi.org/10.1046/j.1365-2869.2003.00337.x\n\n\nDanisch, S., & Krumbiegel, J. (2021). Makie.jl: Flexible high-performance data visualization for julia. Journal of Open Source Software, 6(65), 3349. https://doi.org/10.21105/joss.03349\n\n\nSarkar, D. (2008). Lattice: Mutivariate data visualization with r. Springer-Verlag GmbH. https://www.ebook.de/de/product/11429038/deepayan_sarkar_lattice.html\n\n\n\n\n\n\n Back to top", "crumbs": [ "Worked examples", "Analysis of the sleepstudy data" diff --git a/sitemap.xml b/sitemap.xml index 634b132..c490a45 100644 --- a/sitemap.xml +++ b/sitemap.xml @@ -26,7 +26,7 @@ https://RePsychLing.github.io/SMLP2024/arrow.html - 2024-09-13T10:21:32.364Z + 2024-09-13T11:19:01.984Z https://RePsychLing.github.io/SMLP2024/kb07.html @@ -78,11 +78,11 @@ https://RePsychLing.github.io/SMLP2024/contrasts_kwdyz11.html - 2024-09-13T10:21:36.964Z + 2024-09-13T11:19:01.985Z https://RePsychLing.github.io/SMLP2024/sleepstudy.html - 2024-09-13T10:21:32.365Z + 2024-09-13T11:19:01.985Z https://RePsychLing.github.io/SMLP2024/bootstrap.html diff --git a/sleepstudy.html b/sleepstudy.html index fa95cf3..b634c1b 100644 --- a/sleepstudy.html +++ b/sleepstudy.html @@ -636,7 +636,7 @@

1 Loading the dat
Code -
let f = Figure(; size=(800,400))
+
let f = Figure(; size=(700, 400))
   yrange = maximum(sleepstudy.reaction) - minimum(sleepstudy.reaction)
   xrange = maximum(sleepstudy.days) - minimum(sleepstudy.days)
   
@@ -693,7 +693,7 @@ 
1 Loading the dat
 

 

 

-
+
 

 

 Figure 1: Average response time versus days of sleep deprivation by subject