{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# LMM predictions and prediction intervals\n", "\n", "Below we will fit a linear mixed model using the Ruby gem [mixed\\_models](https://github.com/agisga/mixed_models) and demonstrate the available prediction methods." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data and linear mixed model\n", "\n", "We use the same data and model formulation as in several previous examples, where we have looked at various parameter estimates ([1](http://nbviewer.ipython.org/github/agisga/mixed_models/blob/master/notebooks/LMM_model_fitting.ipynb)) and demostrated many types hypotheses tests as well as confidence intervals ([2](http://nbviewer.ipython.org/github/agisga/mixed_models/blob/master/notebooks/LMM_tests_and_intervals.ipynb)).\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 of *Species* $spcs$ who is at the *Location* $lctn$ as:\n", "\n", "$$Aggression = \\beta_{0} + \\gamma_{spcs} + Age \\cdot \\beta_{1} + b_{lctn,0} + Age \\cdot b_{lctn,1} + \\epsilon,$$\n", "\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*).\n", "\n", "We fit this model in `mixed_models` using a syntax familiar from the `R` package `lme4`." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47316264430760 rows: 5 cols: 4
coefsdz_scoreWaldZ_p_value
intercept1016.286720702345960.1972749576905416.8826034304150770.0
Age-0.065316153427889070.0898848636725299-0.72666465475043740.46743141066211646
Species_lvl_Human-499.693695290208550.2682523406941929-1862.7747813759370.0
Species_lvl_Ood-899.56932135357650.28144708140043684-3196.22899224060030.0
Species_lvl_WeepingAngel-199.588958042007020.27578357795259995-723.71589172837250.0
" ], "text/plain": [ "\n", "#\n", " coef sd z_score WaldZ_p_va \n", " intercept 1016.28672 60.1972749 16.8826034 0.0 \n", " Age -0.0653161 0.08988486 -0.7266646 0.46743141 \n", "Species_lv -499.69369 0.26825234 -1862.7747 0.0 \n", "Species_lv -899.56932 0.28144708 -3196.2289 0.0 \n", "Species_lv -199.58895 0.27578357 -723.71589 0.0 \n" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "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", "model_fit.fix_ef_summary" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Predictions and prediction intervals\n", "\n", "Often, the objective of a statistical model is the prediction of future observations based on new data input.\n", "\n", "We consider the following new data set containing age, geographic location and species for ten individuals." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47316263806300 rows: 10 cols: 3
AgeSpeciesLocation
0209DalekOodSphere
190OodEarth
2173OodAsylum
3153HumanAsylum
4255WeepingAngelOodSphere
5256WeepingAngelAsylum
637DalekEarth
7146WeepingAngelEarth
8127WeepingAngelAsylum
941OodAsylum
" ], "text/plain": [ "\n", "#\n", " Age Species Location \n", " 0 209 Dalek OodSphere \n", " 1 90 Ood Earth \n", " 2 173 Ood Asylum \n", " 3 153 Human Asylum \n", " 4 255 WeepingAng OodSphere \n", " 5 256 WeepingAng Asylum \n", " 6 37 Dalek Earth \n", " 7 146 WeepingAng Earth \n", " 8 127 WeepingAng Asylum \n", " 9 41 Ood Asylum \n" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "newdata = Daru::DataFrame.from_csv '../examples/data/alien_species_newdata.csv'\n", "newdata.vectors = Daru::Index.new(newdata.vectors.map { |v| v.to_sym })\n", "newdata" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Point estimates\n", "\n", "Based on the fitted linear mixed model we can predict the aggression levels for the inidividuals, where we can specify whether the random effects estimates should be included in the calculations or not." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Predictions of aggression levels on a new data set:\n" ] }, { "data": { "text/plain": [ "[1070.9125752531208, 182.45206492790737, -17.06446875476354, 384.7881586199103, 876.1240725686446, 674.7113391148862, 1092.6985606350866, 871.1508855262363, 687.4629975728096, -4.016260100144294]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "puts \"Predictions of aggression levels on a new data set:\"\n", "pred = model_fit.predict(newdata: newdata, with_ran_ef: true)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can add the computed predictions to the data set, in order to see better which of the individuals are likely to be particularly dangerous." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47316262633840 rows: 10 cols: 4
AgeSpeciesLocationPredicted_Agression
0209DalekOodSphere1070.9125752531208
190OodEarth182.45206492790737
2173OodAsylum-17.06446875476354
3153HumanAsylum384.7881586199103
4255WeepingAngelOodSphere876.1240725686446
5256WeepingAngelAsylum674.7113391148862
637DalekEarth1092.6985606350866
7146WeepingAngelEarth871.1508855262363
8127WeepingAngelAsylum687.4629975728096
941OodAsylum-4.016260100144294
" ], "text/plain": [ "\n", "#\n", " Age Species Location Predicted_ \n", " 0 209 Dalek OodSphere 1070.91257 \n", " 1 90 Ood Earth 182.452064 \n", " 2 173 Ood Asylum -17.064468 \n", " 3 153 Human Asylum 384.788158 \n", " 4 255 WeepingAng OodSphere 876.124072 \n", " 5 256 WeepingAng Asylum 674.711339 \n", " 6 37 Dalek Earth 1092.69856 \n", " 7 146 WeepingAng Earth 871.150885 \n", " 8 127 WeepingAng Asylum 687.462997 \n", " 9 41 Ood Asylum -4.0162601 \n" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "newdata = Daru::DataFrame.from_csv '../examples/data/alien_species_newdata.csv'\n", "newdata.vectors = Daru::Index.new(newdata.vectors.map { |v| v.to_sym })\n", "newdata[:Predicted_Agression] = pred\n", "newdata" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Interval estimates\n", "\n", "Since the estimated fixed and random effects coefficients most likely are not exactly the true values, we probably should look at interval estimates of the predictions, rather than the point estimates computed above.\n", "\n", "Two types of such interval estimates are currently available in `LMM`. On the one hand, a *confidence interval* is an interval estimate of the mean value of the response for given covariates (i.e. a population parameter); on the other hand, a *prediction interval* is an interval estimate of a future observation (for further explanation of this distinction see for example )." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "88% confidence intervals for the predictions:\n" ] }, { "data": { "text/html": [ "
Daru::DataFrame:47316259596660 rows: 10 cols: 3
predlower88upper88
01002.6356446359171906.2754736170911098.995815654743
1110.8389455402593717.15393113018095204.5239599503378
2105.4177048057446210.164687937713381200.67072167377586
3506.59965393767027411.8519191795299601.3473886958107
4800.0421435362272701.9091174988788898.1751695735755
5799.9768273827992701.8009453018722898.1527094637263
61013.870023025514920.4439313191591107.296114731869
7807.1616042598671712.571759209002901.7514493107321
8808.402611174997714.191640124036902.613582225958
9114.0394370582259920.614034870631627207.46483924582034
" ], "text/plain": [ "\n", "#\n", " pred lower88 upper88 \n", " 0 1002.63564 906.275473 1098.99581 \n", " 1 110.838945 17.1539311 204.523959 \n", " 2 105.417704 10.1646879 200.670721 \n", " 3 506.599653 411.851919 601.347388 \n", " 4 800.042143 701.909117 898.175169 \n", " 5 799.976827 701.800945 898.152709 \n", " 6 1013.87002 920.443931 1107.29611 \n", " 7 807.161604 712.571759 901.751449 \n", " 8 808.402611 714.191640 902.613582 \n", " 9 114.039437 20.6140348 207.464839 \n" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "puts \"88% confidence intervals for the predictions:\"\n", "ci = model_fit.predict_with_intervals(newdata: newdata, level: 0.88, type: :confidence)\n", "Daru::DataFrame.new(ci, order: [:pred, :lower88, :upper88])" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "88% prediction intervals for the predictions:\n" ] }, { "data": { "text/html": [ "
Daru::DataFrame:47316258683700 rows: 10 cols: 3
predlower88upper88
01002.6356446359171809.91005014591041195.3612391259237
1110.83894554025937-76.53615884686141298.2140499273802
2105.41770480574462-85.09352864481423295.92893825630347
3506.59965393767027317.0988995529618696.1004083223787
4800.0421435362272603.7713980881146996.3128889843398
5799.9768273827992603.6203777073699996.3332770582285
61013.870023025514827.01272323178051200.7273228192475
7807.1616042598671617.9767304115936996.3464781081406
8808.402611174997619.9754792487822996.8297431012118
9114.03943705822599-72.8161447158925300.8950188323445
" ], "text/plain": [ "\n", "#\n", " pred lower88 upper88 \n", " 0 1002.63564 809.910050 1195.36123 \n", " 1 110.838945 -76.536158 298.214049 \n", " 2 105.417704 -85.093528 295.928938 \n", " 3 506.599653 317.098899 696.100408 \n", " 4 800.042143 603.771398 996.312888 \n", " 5 799.976827 603.620377 996.333277 \n", " 6 1013.87002 827.012723 1200.72732 \n", " 7 807.161604 617.976730 996.346478 \n", " 8 808.402611 619.975479 996.829743 \n", " 9 114.039437 -72.816144 300.895018 \n" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "puts \"88% prediction intervals for the predictions:\"\n", "pi = model_fit.predict_with_intervals(newdata: newdata, level: 0.88, type: :prediction)\n", "Daru::DataFrame.new(pi, order: [:pred, :lower88, :upper88])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Remark**: You might notice that `#predict` with `with_ran_ef: true` produces some values outside of the confidence intervals, because the confidence intervals are computed from `#predict` with `with_ran_ef: false`.\n", "However, `#predict` with `with_ran_ef: false` should always give values which lie in the center of the confidence or prediction intervals." ] } ], "metadata": { "anaconda-cloud": {}, "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 }