{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# LMM hypotheses tests and confidence intervals\n", "\n", "Often statistical models are used in order to determine which of the predictor variables have a significant relationship with the response variable. `LMM` has a number of methods to aid with this kind of statistical inference.\n", "\n", "Below we will fit a linear mixed model using the Ruby gem [mixed\\_models](https://github.com/agisga/mixed_models), and demostrate various types of hypotheses tests and confidence intervals available for objects of class `LMM`.\n", "\n", "\n", "## Contents\n", "\n", "* Data and linear mixed model\n", "\n", "* Likelihood ratio test\n", " - Chi squared p-value\n", " - Bootstrap p-value\n", "\n", "* Fixed effects hypotheses tests\n", " - Likelihood ratio test\n", " - Bootstrap test\n", " - Variances and covariances\n", " - Wald Z test\n", "\n", "* Fixed effects confidence intervals\n", " - Bootstrap confidence intervals\n", " - Wald Z confidence intervals\n", " - All types of intervals at once\n", "\n", "* Random effects hypotheses tests\n", " - Likelihood ratio test\n", " - Bootstrap test" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data and linear mixed model\n", "\n", "We use the same data and model formulation as in a [previous example](http://nbviewer.ipython.org/github/agisga/mixed_models/blob/master/notebooks/LMM_model_fitting.ipynb), where we have looked at various parameter estimates, and have shown that the model fit is good.\n", "\n", "The data set, which is simulated, contains two numeric variables *Age* and *Aggression*, and two categorical variables *Location* and *Species*. These data are available for 100 (human and alien) individuals.\n", "\n", "We model the *Aggression* level of an individual as a linear function of the *Age* (*Aggression* decreases with *Age*), with a different constant added for each *Species* (i.e. each species has a different base level of aggression). Moreover, we assume that there is a random fluctuation in *Aggression* due to the *Location* that an individual is at. Additionally, there is a random fluctuation in how *Age* affects *Aggression* at each different *Location*. \n", "\n", "Thus, the *Aggression* level of an individual of *Species* $spcs$ who is at the *Location* $lctn$ can be expressed as:\n", "$$Aggression = \\beta_{0} + \\gamma_{spcs} + Age \\cdot \\beta_{1} + b_{lctn,0} + Age \\cdot b_{lctn,1} + \\epsilon,$$\n", "where $\\epsilon$ is a random residual, and the random vector $(b_{lctn,0}, b_{lctn,1})^T$ follows a multivariate normal distribution (the same distribution but different realizations of the random vector for each *Location*). That is, we have a linear mixed model with fixed effects $\\beta_{0}, \\beta_{1}, \\gamma_{Dalek}, \\gamma_{Ood}, \\dots$, and random effects $b_{Asylum,0}, b_{Asylum,1}, b_{Earth,0},\\dots$.\n", "\n", "We fit the model with the convenient method `LMM#from_formula`, which mimics the behaviour of the function `lmer` from the `R` package `lme4`." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fixed effects terms estimates and some diagnostics:\n", "\n", "#\n", " coef sd z_score WaldZ_p_value \n", " intercept 1016.2867207023459 60.19727495769054 16.882603430415077 0.0 \n", " Age -0.06531615342788907 0.0898848636725299 -0.7266646547504374 0.46743141066211646 \n", " Species_lvl_Human -499.69369529020855 0.2682523406941929 -1862.774781375937 0.0 \n", " Species_lvl_Ood -899.5693213535765 0.28144708140043684 -3196.2289922406003 0.0 \n", "Species_lvl_WeepingA -199.58895804200702 0.27578357795259995 -723.7158917283725 0.0 \n", "\n", "Random effects correlation structure:\n", "\n", "#\n", " Location Location_Age \n", " Location 104.26376362 -0.059863903 \n", "Location_Age -0.059863903 0.1556707755 \n", "\n" ] } ], "source": [ "require 'mixed_models'\n", "\n", "alien_species = Daru::DataFrame.from_csv '../examples/data/alien_species.csv'\n", "# mixed_models expects that all variable names in the data frame are ruby Symbols:\n", "alien_species.vectors = Daru::Index.new(alien_species.vectors.map { |v| v.to_sym })\n", "\n", "model_fit = LMM.from_formula(formula: \"Aggression ~ Age + Species + (Age | Location)\", \n", " data: alien_species)\n", "\n", "puts \"Fixed effects terms estimates and some diagnostics:\"\n", "puts model_fit.fix_ef_summary.inspect(20)\n", "puts \"Random effects correlation structure:\"\n", "puts model_fit.ran_ef_summary.inspect(12)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Likelihood ratio test\n", "\n", "Given two nested models, the `LMM.likelihood_ratio_test` tests whether the restricted simpler model is adequate. In this context 'nested' means that all predictors used in the restricted model must also be predictors in the full model (i.e. one model is a reduced version of the other more complex model). This method works only if both models were fit using the deviance (as opposed to REML criterion) as the objective function for the minimization (i.e. fit with `reml: false`). `LMM.likelihood_ratio_test` returns the p-value of the test.\n", "\n", "Two methods are available to compute the p-value: approximation by a Chi squared distribution as delineated in section 2.4.1 in Pinheiro & Bates \"Mixed Effects Models in S and S-PLUS\" (2000), and a method based on bootstrapping as delineated in section 4.2.3 in Davison & Hinkley \"Bootstrap Methods and their Application\" (1997).\n", "\n", "For example we can test the model formulation as above against a simpler model, which assumes that *Age* is neither a fixed effect nor a random effect. We can compute a likelihood ratio using the method `LMM.likelihood_ratio`." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "454.3606613877624" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "complex_model = LMM.from_formula(formula: \"Aggression ~ Age + Species + (Age | Location)\", \n", " data: alien_species, reml: false)\n", "simple_model = LMM.from_formula(formula: \"Aggression ~ Species + (1 | Location)\", \n", " data: alien_species, reml: false)\n", "\n", "LMM.likelihood_ratio(simple_model, complex_model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Chi squared p-value\n", "\n", "We perform the likelihood ratio test using the method `LMM.likelihood_ratio_test` with `method: :chi2` to use the Chi squared approximation for the p-values." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3.693825760412622e-98" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "chi2_p_value = LMM.likelihood_ratio_test(simple_model, complex_model, method: :chi2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The p-value is tiny, which implies that *Age* is a significant predictor of *Aggression*, and that the more complex model should be prefered." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bootstrap p-value\n", "\n", "However, often one may not be sure whether the assumptions required for the validity of the Chi squared test are satisfied. In that case, one can compute a p-value with the bootstrap method by specifying `method: :bootstrap`, which by default uses all available CPUs in parallel." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.000999000999000999" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bootstrap_p_value = LMM.likelihood_ratio_test(simple_model, complex_model, method: :bootstrap, nsim: 1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Even though the p-value is not as extreme as with the Chi squared test, it still shows significance of the variable *Age*." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fixed effects hypotheses tests\n", "\n", "Significance tests for the fixed effects can be performed with `LMM#fix_ef_p` (or its alias `LMM#fix_ef_test`). For a given fixed effects coefficient estimate the tested null hypothesis is that the true value of the coefficient is zero (i.e. no linear relationship to the response). \n", "\n", "That is, for the above model formulation we carry out hypotheses tests for each fixed effects terms $\\beta_{i}$ or $\\gamma_{species}$, testing the null $H_{0} : \\beta_{i} = 0$ against the alternative $H_{a} : \\beta_{i} \\neq 0$, or respectively the null $H_{0} : \\gamma_{species} = 0$ against the alternative $H_{a} : \\gamma_{species} \\neq 0$.\n", "\n", "`LMM` currently offers three methods of hypotheses testing for fixed effects: *Wald Z test*, *likelihood ratio test*, and a *bootstrap test*. For a good discussion of the validity of different methods see [this entry from the wiki of the r-sig-mixed-models mailing list](http://glmm.wikidot.com/faq). Moreover, due to the equivalence of hypotheses tests and confidence intervals, an additional hypothesis testing tool are bootstrap confidence intervals which are described in a different section below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Likelihood ratio test for the fixed effects\n", "\n", "The likelihood ratio test for fixed effects is actually merely a convenience method. It is a convenient interface to `LMM.likelihood_ratio_test` with `method: :chi2` described above. For example, we can test whether the fixed effects term *Age* is a significant predictor of *Aggression* in `complex_model` as follows." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.40176699624310697" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lrt_p_value = complex_model.fix_ef_p(variable: :Age, method: :lrt)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that *Age* does not seem be significant as a fixed effects term, even though it is in general a significant predictor as we have seen before." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bootstrap test for the fixed effects\n", "\n", "Like the likelihood ratio test, the bootstrap test is merely a convenient shortcut to `LMM.likelihood_ratio_test` with `method: :bootstrap`. We can test the significance of the predictor variable *Age* in `complex_model` with:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.5314685314685315" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bootstrap_p_value = complex_model.fix_ef_p(variable: :Age, method: :bootstrap, nsim: 1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This result confirms the conclusion based of `method: :lrt`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Variances and covariances of the fixed effects coefficient estimates\n", "\n", "The covariance matrix of the fixed effects estimates is returned by `LMM#fix_ef_cov_mat`, and the standard deviations of the fixed effects coefficients are returned by `LMM#fix_ef_sd`." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{:intercept=>60.19727495769054, :Age=>0.0898848636725299, :Species_lvl_Human=>0.2682523406941929, :Species_lvl_Ood=>0.28144708140043684, :Species_lvl_WeepingAngel=>0.27578357795259995}" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model_fit.fix_ef_sd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Wald Z tests for the fixed effects\n", "\n", "Based on the values returned by `LMM#fix_ef_sd`, the [Wald Z test statistics](https://en.wikipedia.org/wiki/Wald_test#Test_on_a_single_parameter) for all fixed effects coefficients are computed by the method `LMM#fix_ef_z`." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{:intercept=>16.882603430415077, :Age=>-0.7266646547504374, :Species_lvl_Human=>-1862.774781375937, :Species_lvl_Ood=>-3196.2289922406003, :Species_lvl_WeepingAngel=>-723.7158917283725}" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model_fit.fix_ef_z" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Based on the above Wald Z test statistics, hypotheses tests for each fixed effects term can be carried out.\n", "The corresponding (approximate) p-values are obtained with `LMM#fix_ef_p` as follows." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{:intercept=>0.0, :Age=>0.46743141066211646, :Species_lvl_Human=>0.0, :Species_lvl_Ood=>0.0, :Species_lvl_WeepingAngel=>0.0}" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model_fit.fix_ef_p(method: :wald)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that *Aggression* of each *Species* is significantly different from the base level (which is the species *Dalek* in this model), while statistically we don't have enough evidence to conclude that *Age* of an individual is a good predictor of his/her/its aggression level (which agrees with the conclusion obtained above with `method: :lrt` and `method: :bootstrap`)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Confidence intervals for the fixed effects\n", "\n", "Confidence intervals for the fixed effects terms can be computed with the method `LMM#fix_ef_conf_int`. The following types of confidence intervals are available:\n", "\n", "* Wald Z intervals\n", "* Basic bootstrap intervals\n", "* Bootstrap normal intervals\n", "* Bootstrap percentile intervals\n", "* Studentized bootstrap intervals" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bootstrap confidence intervals for the fixed effects\n", "\n", "A detailed description of the bootstrap methods available in `LMM#fix_ef_conf_int` is given in [this blog post](http://agisga.github.io/bootstap_confidence_intervals/).\n", "\n", "For example we can compute the studentized bootstrap confidence intervals (also called bootstrap-t) with 98% coverage as follows." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{:intercept=>[874.4165005038727, 1139.1584607682103], :Age=>[-0.2706718136184743, 0.1493216868271593], :Species_lvl_Human=>[-500.28872315642866, -499.09663075641373], :Species_lvl_Ood=>[-900.2572251446181, -898.9341708622952], :Species_lvl_WeepingAngel=>[-200.24313298457852, -198.95081041642916]}" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bootstrap_t_intervals = model_fit.fix_ef_conf_int(level: 0.98, method: :bootstrap, \n", " boottype: :studentized, nsim: 1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that *Age* is the only fixed effects predictor whose confidence interval contains zero, which implies that it probably has little linear relationship as a fixed effect to the response variable *Aggression*." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Wald Z confidence intervals for the fixed effects\n", "\n", "We can use the Wald Z statistic to compute confidence intervals as well. For example 90% confidence intervals for each fixed effects coefficient estimate can be computed as follows." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{:intercept=>[917.2710134723027, 1115.302427932389], :Age=>[-0.21316359921454495, 0.0825312923587668], :Species_lvl_Human=>[-500.1349311310106, -499.2524594494065], :Species_lvl_Ood=>[-900.0322606117453, -899.1063820954076], :Species_lvl_WeepingAngel=>[-200.04258166587707, -199.13533441813698]}" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "conf_int = model_fit.fix_ef_conf_int(level: 0.9, method: :wald)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For greated visual clarity we can put the coefficient estimates and the confidence intervals into a `Daru::DataFrame`:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47402176301980 rows: 5 cols: 3
lower90upper90coef
intercept917.27101347230271115.3024279323891016.2867207023459
Age-0.213163599214544950.0825312923587668-0.06531615342788907
Species_lvl_Human-500.1349311310106-499.2524594494065-499.69369529020855
Species_lvl_Ood-900.0322606117453-899.1063820954076-899.5693213535765
Species_lvl_WeepingAngel-200.04258166587707-199.13533441813698-199.58895804200702
" ], "text/plain": [ "\n", "#\n", " lower90 upper90 coef \n", " intercept 917.271013 1115.30242 1016.28672 \n", " Age -0.2131635 0.08253129 -0.0653161 \n", "Species_lv -500.13493 -499.25245 -499.69369 \n", "Species_lv -900.03226 -899.10638 -899.56932 \n", "Species_lv -200.04258 -199.13533 -199.58895 \n" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = Daru::DataFrame.rows(conf_int.values, order: [:lower90, :upper90], index: model_fit.fix_ef_names)\n", "df[:coef] = model_fit.fix_ef.values\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### All types of confidence intervals at once\n", "\n", "With `method: :all`, `LMM#fix_ef_conf_int` returns a Daru::DataFrame containing the confidence intervals obtained by each of the available methods. The data frame can be printed in form of a nice looking table for inspection. \n", "\n", "For example for the alien species data we obtain all types of 95% confidence intervals with:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47402185048160 rows: 5 cols: 5
interceptAgeSpecies_lvl_HumanSpecies_lvl_OodSpecies_lvl_WeepingAngel
wald_z[898.3022305977157, 1134.2712108069761][-0.2414872478168185, 0.11085494096104034][-500.2194602132623, -499.1679303671548][-900.1209474930289, -899.017695214124][-200.12948391874875, -199.0484321652653]
boot_basic[901.8226468130745, 1142.3725792755527][-0.23613380807522952, 0.12082106157921091][-500.23480713878746, -499.1808531803127][-900.1538925933955, -899.0066844647126][-200.18820176543932, -199.02291203098875]
boot_norm[900.8638839567735, 1134.5288460466913][-0.2473842439597385, 0.11510977545973858][-500.2315699405288, -499.1525612096204][-900.1465247092458, -898.9995195976951][-200.1708564794685, -199.01484568123584]
boot_t[901.8226468130744, 1142.3725792755527][-0.23613380807522952, 0.12082106157921088][-500.2348071387875, -499.1808531803127][-900.1538925933955, -899.0066844647126][-200.18820176543932, -199.02291203098875]
boot_perc[890.2008621291391, 1130.7507945916173][-0.25145336843498906, 0.10550150121945137][-500.2065374001044, -499.15258344162964][-900.1319582424403, -898.9847501137574][-200.1550040530253, -198.98971431857473]
" ], "text/plain": [ "\n", "#\n", " intercept Age Species_lv Species_lv Species_lv \n", " wald_z [898.30223 [-0.241487 [-500.2194 [-900.1209 [-200.1294 \n", "boot_basic [901.82264 [-0.236133 [-500.2348 [-900.1538 [-200.1882 \n", " boot_norm [900.86388 [-0.247384 [-500.2315 [-900.1465 [-200.1708 \n", " boot_t [901.82264 [-0.236133 [-500.2348 [-900.1538 [-200.1882 \n", " boot_perc [890.20086 [-0.251453 [-500.2065 [-900.1319 [-200.1550 \n" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ci = model_fit.fix_ef_conf_int(method: :all, nsim: 1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since here we are dealing with data that was simulated according to the assumptions of the linear mixed model, all parameters end up approximately meeting the normality assumptions, and therefore all confidence interval methods turn out to be pretty much equivalent. Often when analyzing less ideal data, this will not be the case. Then it might be necessary to compare different types of confidence intervals in order to draw the right conclusions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Random effects hypotheses tests\n", "\n", "We can test individual random effects for siginificance with the method `LMM#ran_ef_p` (or its alias `LMM#ran_ef_test`). It offers two methods, `:lrt` and `:bootstrap`. Both are in fact merely convenient interfaces to `LMM.likelihood_ratio_test` described above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Likelihood ratio test for individual random effects\n", "\n", "The likelihood ratio test for random effects can only be performed if the model was fit with option `reml: false`. \n", "\n", "For example we can test the random intercept term (due to *Location*) of `complex_model` as follows." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1.3846568414031767e-151" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "complex_model.ran_ef_p(variable: :intercept, grouping: :Location, method: :lrt)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This suggests that the variability in *Aggression* is significantly influenced by the effect of *Location*." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bootstrap test for individual random effects\n", "\n", "This test can also be performed only if the model was fit with`reml: false`.\n", "\n", "We can test whether there is a significant random variation of the effect of *Age* due to *Location* using `LMM#ran_ef_p` with `method: :bootstrap` as follows." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.000999000999000999" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "complex_model.ran_ef_p(variable: :Age, grouping: :Location, method: :bootstrap, nsim: 1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that in fact the change of *Age* due to *Location* is a significant random effect." ] } ], "metadata": { "kernelspec": { "display_name": "Ruby 2.3.1", "language": "ruby", "name": "ruby" }, "language_info": { "file_extension": ".rb", "mimetype": "application/x-ruby", "name": "ruby", "version": "2.3.1" } }, "nbformat": 4, "nbformat_minor": 0 }