{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Data Analysis in Ruby\n", "\n", "## Introduction\n", "\n", "This is an example of a data analysis in Ruby using [`daru`](https://github.com/v0dro/daru) to organize and manipulate data, [`mixed_models`](https://github.com/agisga/mixed_models) to fit a statistical model, and [`gnuplotrb`](https://github.com/dilcom/gnuplotrb) for visualization.\n", "\n", "We will analyze [these data](http://archive.ics.uci.edu/ml/datasets/BlogFeedback)1 from the UCI machine learning repository2, which originate from blog posts from various sources in 2010-2012.\n", "\n", "Here, we use the data to investigate the following question:\n", "\n", "> *If a blog post has received a comment in the first 24 hours after publication, how many more comments will it receive before 24 hours after its publication have passed?*\n", "\n", "After doing the entire data cleaning and preprocessing with `daru`, we use `mixed_models` to fit a linear mixed model. The fitted model can be used to make predictions for future observations, and make inferences about the relationships between different variables in the data.\n", "\n", "------------------------------------------------------------------------\n", "\n", "[1] For more information on the data set we refer to: Buza, K. (2014). *Feedback Prediction for Blogs*. In Data Analysis, Machine Learning and Knowledge Discovery (pp. 145-152). Springer International Publishing.\n", "\n", "[2] Lichman, M. (2013). *UCI Machine Learning Repository* . Irvine, CA: University of California, School of Information and Computer Science." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data preprocessing with `daru`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since `daru` requires CSV files to have a header line, we add a header to the data file and save the new data frame." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "without_header = '../examples/data/blogData_train.csv'\n", "with_header = '../examples/data/blogData_train_with_header.csv'\n", "colnames = (1..281).to_a.map { |x| \"v#{x}\" }\n", "header = colnames.join(',')\n", "File.open(with_header, 'w') do |fo|\n", " fo.puts header\n", " File.foreach(without_header) do |li|\n", " fo.puts li\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can load the data with `daru`." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of rows: 52397\n", "Number of columns: 281\n" ] } ], "source": [ "require 'daru'\n", "df = Daru::DataFrame.from_csv '../examples/data/blogData_train_with_header.csv'\n", "puts \"Number of rows: #{df.nrows}\"\n", "puts \"Number of columns: #{df.ncols}\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Selecting and renaming data columns\n", "\n", "We don't want to keep most of the variables for further analysis, as the great majority of the columns represent bag-of-words features and summary statistics for other attributes.\n", "So, we select the data columns which we want to keep, and assign them meaningful names." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47139036762220 rows: 10 cols: 13
host_comments_avghost_trackbacks_avgcommentslengthmotuwethfrsasuparentsparents_comments
034.5675660.9729732.00.00.00.00.01.00.00.00.00.00.0
134.5675660.9729735.00.00.00.01.00.00.00.00.00.00.0
234.5675660.9729735.00.00.00.01.00.00.00.00.00.00.0
334.5675660.9729732.00.00.00.00.01.00.00.00.00.00.0
434.5675660.9729732.00.00.00.00.01.00.00.00.00.00.0
534.5675660.9729735.00.00.00.01.00.00.00.00.00.00.0
634.5675660.9729735.00.00.00.01.00.00.00.00.00.00.0
734.5675660.9729732.00.00.00.00.01.00.00.00.00.00.0
834.5675660.9729732.00.00.00.00.01.00.00.00.00.00.0
934.5675660.9729732.00.00.00.00.01.00.00.00.00.00.0
" ], "text/plain": [ "\n", "#\n", " host_comme host_track comments length mo tu we th fr sa su parents parents_co \n", " 0 34.567566 0.972973 2.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 1 34.567566 0.972973 5.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 2 34.567566 0.972973 5.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 3 34.567566 0.972973 2.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 4 34.567566 0.972973 2.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 5 34.567566 0.972973 5.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 6 34.567566 0.972973 5.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 7 34.567566 0.972973 2.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 8 34.567566 0.972973 2.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 9 34.567566 0.972973 2.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# select a subset of columns of the data frame\n", "keep = ['v16', 'v41', 'v54', 'v62', 'v270', 'v271', 'v272', \n", " 'v273', 'v274', 'v275', 'v276', 'v277', 'v280']\n", "blog_data = df[*keep]\n", "df = nil\n", "\n", "# assign meaningful names for the selected columns\n", "meaningful_names = [:host_comments_avg, :host_trackbacks_avg, \n", " :comments, :length, :mo, :tu, :we, :th, \n", " :fr, :sa, :su, :parents, :parents_comments]\n", "blog_data.vectors = Daru::Index.new(meaningful_names)\n", "\n", "# the resulting data set\n", "blog_data.head" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Selecting a subset of the data rows\n", "\n", "As can be seen in the above output, the length of the text in a blog post is often given as zero. Those are probably missing values, and we get rid of the corresponding observations. \n", "\n", "We also delete observation which have zero comments in the first 24 hours after publication, to comply with our research objective stated in the beginning." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Remaining number of rows: 22435\n" ] } ], "source": [ "nonzero_ind = blog_data[:length].each_index.select do |i| \n", " blog_data[:length][i] > 0 && blog_data[:comments][i] > 0\n", "end\n", "blog_data = blog_data.row[*nonzero_ind]\n", "puts \"Remaining number of rows: #{blog_data.nrows}\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Replacing and transforming variables\n", "\n", "##### Creating a categorical \"day\" variable\n", "\n", "For a more clear representation of the data, and in order to use the day of the week as a grouping variable for the observations, we replace the respective seven 0-1-valued columns with one column of categorical data (with values 'mo', 'tu', 'we', 'th', 'fr', 'sa', 'su')." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47139045191780 rows: 3 cols: 7
host_comments_avghost_trackbacks_avgcommentslengthparentsparents_commentsday
1221110.300870.074.03501.00.00.0we
1222110.300870.074.03501.00.00.0we
1223110.300870.0218.04324.00.00.0th
" ], "text/plain": [ "\n", "#\n", " host_comme host_track comments length parents parents_co day \n", " 1221 110.30087 0.0 74.0 3501.0 0.0 0.0 we \n", " 1222 110.30087 0.0 74.0 3501.0 0.0 0.0 we \n", " 1223 110.30087 0.0 218.0 4324.0 0.0 0.0 th \n" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "days = Array.new(blog_data.nrows) { :unknown }\n", "[:mo, :tu, :we, :th, :fr, :sa, :su].each do |d|\n", " ind = blog_data[d].to_a.each_index.select { |i| blog_data[d][i]==1 }\n", " ind.each { |i| days[i] = d.to_s }\n", " blog_data.delete_vector(d)\n", "end\n", "blog_data[:day] = days\n", "blog_data.head 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Replacing two highly correlated variables\n", "\n", "The variable `parents` denotes the number of parent blog posts, where we consider a blog post P as a\n", "parent of blog post B if B is a reply (trackback) to blog post P. \n", "Related to it, the variable `parents_comments` denotes the number of comments that the parents of a blog post received on average.\n", "Clearly, the two variables are highly correlated (as zero `parents` implies zero `parents_comments`). Therefore, we shouldn't include both these variables in the linear mixed model.\n", "\n", "We combine the variables `parents` and `parents_comments` into one variable called `has_parent_with_comments`, which is a categorical variable designating whether a blog post has at least one parent post with at least one comment." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47139031507220 rows: 3 cols: 6
host_comments_avghost_trackbacks_avgcommentslengthdayhas_parent_with_comments
1221110.300870.074.03501.0weno
1222110.300870.074.03501.0weno
1223110.300870.0218.04324.0thno
" ], "text/plain": [ "\n", "#\n", " host_comme host_track comments length day has_parent \n", " 1221 110.30087 0.0 74.0 3501.0 we no \n", " 1222 110.30087 0.0 74.0 3501.0 we no \n", " 1223 110.30087 0.0 218.0 4324.0 th no \n" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# create a binary indicator vector specifying if a blog post has at least \n", "# one parent post which has comments\n", "hpwc = (blog_data[:parents] * blog_data[:parents_comments]).to_a\n", "blog_data[:has_parent_with_comments] = hpwc.map { |t| t == 0 ? 'no' : 'yes'} \n", "blog_data.delete_vector(:parents)\n", "blog_data.delete_vector(:parents_comments)\n", "blog_data.head 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Log transforms\n", "\n", "Some prior experimentation with the data suggests that it is necessary to take $log$-transforms of the response variable `comments` and the predictor variable `host_comments_avg`.\n", "This results in a much better agreement with the normality assumption on the model residuals." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47139031510440 rows: 3 cols: 8
host_comments_avghost_trackbacks_avgcommentslengthdayhas_parent_with_commentslog_commentslog_host_comments_avg
1221110.300870.074.03501.0weno4.304065093204174.7032118138076795
1222110.300870.074.03501.0weno4.304065093204174.7032118138076795
1223110.300870.0218.04324.0thno5.3844950627890894.7032118138076795
" ], "text/plain": [ "\n", "#\n", " host_comme host_track comments length day has_parent log_commen log_host_c \n", " 1221 110.30087 0.0 74.0 3501.0 we no 4.30406509 4.70321181 \n", " 1222 110.30087 0.0 74.0 3501.0 we no 4.30406509 4.70321181 \n", " 1223 110.30087 0.0 218.0 4324.0 th no 5.38449506 4.70321181 \n" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "log_comments = blog_data[:comments].to_a.map { |c| Math::log(c) }\n", "log_host_comments_avg = blog_data[:host_comments_avg].to_a.map { |c| Math::log(c) }\n", "blog_data[:log_comments] = log_comments\n", "blog_data[:log_host_comments_avg] = log_host_comments_avg\n", "blog_data.head 3" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Linear mixed model\n", "\n", "The logarithm of the number of comments of a given blog post is modeled as a function of the blog post text length, the average number of comments and trackbacks per blog article at the hosting website of the blog, and the existence of commented parent blog posts. Additionally, we model random fluctuations of the number of comments due to the day of the week when the blog post was released.\n", "\n", "That is, the number of comments of the $i$th blog post, which is published on weekday $d$, is estimated as\n", "\n", "$$log(comments_i) \\approx \\beta_0 + length_i \\cdot \\beta_1 + log(host\\_comments_i) \\cdot \\beta_2 + host\\_trackbacks_i \\cdot \\beta_3 + parent\\_with\\_comments\\_yes \\cdot \\beta_4 + b_d.$$\n", "\n", "*Note:* It is questionable whether a linear mixed model is a good choice here, because the response variable represents count data, which for obvious reasons cannot follow a normal distribution. However, we can still get a reasonble model fit, and certainly gain some insight into the data from the linear mixed model, even if it wont match the data very well. In fact, by taking the $log$-transform, we partially account for the fact that the response variable represents count data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We fit the model with `mixed_models` and display the estimated fixed effects coefficients (the $\\beta$'s from the above equation)." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{:intercept=>0.18944500002662124, :log_host_comments_avg=>0.7861228519075351, :host_trackbacks_avg=>0.05579925647682859, :length=>2.580329980620255e-05, :has_parent_with_comments_lvl_yes=>-0.48841082096511457}" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "require 'mixed_models'\n", "model_fit = LMM.from_formula(formula: \"log_comments ~ log_host_comments_avg + host_trackbacks_avg + length + has_parent_with_comments + (1 | day)\", data: blog_data)\n", "model_fit.fix_ef" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47139041993940 rows: 5 cols: 4
coefsdz_scoreWaldZ_p_value
intercept0.189445000026621240.019300366951039869.8156164858002540.0
log_host_comments_avg0.78612285190753510.0053227324152766145147.691597205094780.0
host_trackbacks_avg0.055799256476828590.0081028499656306546.8863741416302615.723421736547607e-12
length2.580329980620255e-052.0833561682273268e-0612.3854481531873170.0
has_parent_with_comments_lvl_yes-0.488410820965114570.08884797074497962-5.4971522351028193.859734754030342e-08
" ], "text/plain": [ "\n", "#\n", " coef sd z_score WaldZ_p_va \n", " intercept 0.18944500 0.01930036 9.81561648 0.0 \n", "log_host_c 0.78612285 0.00532273 147.691597 0.0 \n", "host_track 0.05579925 0.00810284 6.88637414 5.72342173 \n", " length 2.58032998 2.08335616 12.3854481 0.0 \n", "has_parent -0.4884108 0.08884797 -5.4971522 3.85973475 \n" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model_fit.fix_ef_summary" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Assess the quality of the fit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Fitted vs. residuals plot\n", "\n", "We can assess the goodness of the model fit (to some extent) by plotting the residuals agains the fitted values. We definately see a pattern here -- the model seems to make somewhat better predictions for observations with a small number of comments. Additionally, a subset of the observations displays what appears to be a linear relationship between the fitted values and the residuals, which is a reason for concern. However, this is due to the fact the the response variable has a discrete range ($log$ transfrom of count data), with possible values at the low end being especially far apart." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAMAAAACDyzWAAABMlBMVEX///8AAACgoKD/AAAAwAAAgP/AAP8A7u7AQADIyABBaeH/wCAAgEDAgP8wYICLAABAgAD/gP9//9SlKir//wBA4NAAAAAaGhozMzNNTU1mZmZ/f3+ZmZmzs7PAwMDMzMzl5eX////wMjKQ7pCt2ObwVfDg///u3YL/tsGv7u7/1wAA/wAAZAAA/38iiyIui1cAAP8AAIsZGXAAAIAAAM2HzusA////AP8AztH/FJP/f1DwgID/RQD6gHLplnrw5oy9t2u4hgv19dyggCD/pQDugu6UANPdoN2QUEBVay+AFACAFBSAQBSAQICAYMCAYP+AgAD/gED/oED/oGD/oHD/wMD//4D//8DNt57w//Cgts3B/8HNwLB8/0Cg/yC+vr5fX18fHx8/Pz/f39+/v7+fn595gRnyAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4nO1dC3bjKgxV19H9sI6kafe/hZk0BvQFYeM4TnXfeZ0EhABxLQnsJACBQCAQCAQCgUAgEAgEAoFAIBAIBAKBQCAQCAQCgcBfxOXjjuvRwwj8VXzdjh5B4E8jCBg4FNcgYOBIXL+/Pq4/R48i8Gfx9XmD72Bg4FB8s13wRyDQxGwGfjTfbtL1tLbR8Tk6Zvjdg3x/TdR/SqtEx89rTKHmgH/PKtHx8xpT3C5XuQv+e1aJjp/X+AX0B06OIGDgUAQBA4ciCBg4FEHAwKEIAgYORRAwcCiCgIFDEQQMHIogYOBQBAEDhyIIGDgUQcDAoQgCBg5FEDBwKIKAgUMRBAwciiBg4FAEAQOHIggYOBRBwMChCAIGDkUQMHAogoCBQxEEDByKIGDgUAQBA4ciCBg4FEHAwKGYTJAfri8IGGhiLkFu1yBgYAhzCfL1HQQMDGEqQT5/hL4gIEM6egAvhpkEuX8/dBCwiXTH0YN4Kcz8it77TyUFAZtIED6QYiJBLp+Kvl1+CuK0+PV+4QIzJpNj+eGRT1o4T/87IDwgx2yCRAhuInJAjiDgkxH0owgCBg5F3AueAOnVws95EQTcDJnXRabnRxBwM9Lvf6Qk7eoD34rbQcCtSIk7PFkyvb+ddB+AIOBmiLO9Xwe4G0fe7CgxCLgVwt/t7QDhrW6mBAE3Q+SAuzrA8ICvpf8giJQPeMF+IThywFfSfwjErkPUxy7YiyDgCnSi4LulabsiCDiOLsGenKadmupBwBXoEWw0TdvEoJPnhEHAFeiv+RD9tjHo5LviIOAqzFzxbQw6e8IZBDwabQb12Rke8FD9Z0Qqf8pbg0Ep1dsqJkUjBzxS/ytgbPkTZlUpsDU//m/R7MT0CwJux5otL3FptoLHUw35z7mJZiEIuBUaNVyn1Jl1rRCc/z/7VsNGEHAjFGq0fSL3gDa3EorW4QFfVP/xkNRokgWzqiOegPjAt3SAQcDNUB6GAfN5hMq8RMtUxYBd41vSLwg4A1oGSPetiVb1FAD1fm+NIOBmCF+3ZG2lvFLRu5fIGpph9z24GQTcCJbQLYWlfHlXdhzgIo7D/3UPBk/CzyDgRiQj30O3MJY//r2Ex1Ha/BS7nJdGEHAblkjJ6UI4UJjnTuq6gq2Dwd/PqJwmf3wPAh5oas0DouO7cfIVBZ3q2qnKfTjJwfU7EPDAcKOHO3p+DFBckn+UjdMXdpIoe889nYF/Uwnyef24/uyo38KBxl665j4IyCYES47qV66tskFOqmLGzynY0bpTv6L3BrcvxsAnEPDAcKN3rZ5Er2EECuNEe8qPKChDSPlzyhPpt2eEmUiQ7/uf23U3/SaO94C8FIdfXDSo3AjwWZniAXfZAO9q39kEOYSAh+eAein5vg4YW8QSZRsZJj/mhnoAOTqNZuWuEWYyQW6X54dgOC7dTqYPZOXaMzMtvYuLK3uZJOrYI9W1GyVq92bRoeyJPODHx2VX/a+F4oaU5ZOp2UhKmAqXMhNJnfL0TXlZvObAPIRGVn+SHPAO6QHf+HdCmucdrJwvIq4WzesuRonladmC5KY0BKdxAjpC7E7024McR+SAB4EeMtMq43iOtwXQ/UvJ/VQC0vQP16P3g2eORyUxMwny6/v+DgEXZ6M6MJUCaNO6CNQCLEEdH++AEBCY/8LtRlzgjiG2h+edAx42xZ1gR7qk1eLgiKMkPswjng3QYR9UjnKOpXLrlznOsbkchpke6nL9uH4b+g+9yvaA7WbqSTHdHpS9MSVJkasv8AMOhJKLjsS9KWfniYz9rHvBe+UZx9nZnpEWB7OPy39ToVF1d4VziXJJy/vw68o23Ganac/Gkwi402HmkZe63XdlQEJloHjA6sioZyM0Wo5WUDBOtKuFwUz1SRh4bg946P7N6jmHYMUzlYCZaITmjCTva0ul/4Sq0fvT8O95BNzDJm6/uvdqJPGOZHaAxpk4v0oah7fUlKbqfidR1D6Fc509v6l43vOAKyZhNUn4Rd8L7u0QuH70Prsn9bhY7GwTaVmppPo0Gp1pAK4Edk+8Iber/V74gVRr3sr6tg20d6CW+vkVghO8uiGphycLj8iZ4CJQ8kDeKaEcLKKUyP57IkVONbhRPgWvTECw7IEp2BDM8vtsgKr6TACt6tE1lUE74vwWCTBmKRdija6Jk447U3zbrjkJqNFbzgJ2tN/rEtBybcKsPQPxE7nJoN4NdVUIlAtZCMZvkP9BSWMhEemven38iiWXWQHgDyjbc8gXCmiif9cDyt+h/C0nhxLQMRBe3y1jsfVrfCrDKm+ETFVdh1ipA9V/Kv5PnucIT1iH4gnEi631S3mC/Wy8MgEXuxrlvMDSAtvjR78D8hYN3GZdLc41OAdD5yrs+uHcLH2l2ndu6k2SAQ/BygJ3wgsTsGEPzSsYKqzMZsY4lMEs6w1yNQsl2JjlqWFRmpSEg3hA5gU5AHXZtkA+BNrT1+k4koDdqdr28FtpWYAtZm0vIM0S8soLv8PcWa2kfKlq8P4CKMVLwCUej3GPOtk2AVHcd1lkIo4joOtq226PSogNOuyRKMxZll8wEJbQWRY6kTbCAVaeJcqNpUpxf7ikdEoGaU3DrtoZBxIQBue81kDbr+zGtcLjLKEFqcOeijpVzQHWRBJTi/XZJGBieWRzIo0Z7ozDCOhITJj8ahNNsK4Zf/kscoampV45G+AeUIwQc7WSqfozqBRVKbj8QWln3RO1rqTGRHfDC3hA35w3sWg/s/JcrzBHbtZLdcIvtBEiKUa1IptKTDcAZasMWBXyu+psnu4ID88BnXMedZjTYWyGst+pMvklYibNx0hLhQ2EQsgBJtK0ujri/DABQVaUAehml77ZYYNtOHwX7EzQugTcl5vaei3RltKFNkBxj2qgTUQd5HLSGmk3XZ+WF7J8Eg1eeF4A3lWV0Tm7FUefA7o9W5n+IUm0vEy4t9H3I7WtkvRlftC6TBmcCuZAjAfEWCecITA1aLSMamhY1lB1G8zA0QR0Tyuvhc60zcZpN1Yuk7JWdc2JTFnsvNa4FhEDsBj2qjUD5Dlg0c45xy6Kwi/UOvHhl86QWjqR/P8+SdDxBCx27c2tcRluNY5Ba9q3Ga3Q2nIfV9syL4c8E17jwiTmYsUjBZX0tR9GQO2kkA0v8TkwV56wad/UA5a59jiwSMHiFaSODcbpNpejy7TCTgPJpLLrtHNASgJgq89dl94+D4XrJBG8xG88hUJR4JW0yzog1wqN4gUI+ICTQtXx0OJNxuk5UE07plWuRXQgnJMenvukJIZBg2r2ocQ18TElgoXCtUKZY50ErqTXA6Hj+xLQG0SJhUmF3aQv1mG/EfeXAZHKuvbVf2iqi6dSnBLTJFigWUvsQqowY6Y6Du1OS+mEcHo2XoWA/iBa9nZOvdjo1kXcvrjbFwdtS70QO/ug0+BpHWZMeYNog6iIyMCTUcIzGZyVDkWjHJRpVjBkdT9eh4DY+F1Zf8zV76IafZs6vA4Sp+xlI6C0NhmBKgHdCl7UCroUmqtUwt2BVlWoxhiYamklpBjsFLwMAYnxe5LZ7A6l5VQNdFeGL/Nmf+3B5KsnoT9olLx94aamujYgLEVEyHOXxFMpWMt4eVWrELA4WchUfeMcMCO5KGhd6booFK4KplUtdogtfwz9QJas6OMeBHfAPJqpsUewlKq9CJGKUkwphZl1xEwrtseAtcfxYgS0FkUKKimUoTHV0CE0L8ESTE2dsZTFyyuGFwsnbawHvPKC+zTZQ/wCFki7DC1ta2GSQ6h1VqbSCxEbMPt3Qqxvx/ICXY4OOWcURspYJESr4GjM62guhxP/BJSb3MfSdppTVhml8avaS+aIiYyyTAa55MSa0u3zVjiUvNrvhCBLdcMweC/LRgBF/9udqLUaQajqmn6WZSfV9TkBPJrCA7K5qA4/lfM9xkQxJuCkZPlCJj8RSAN27cDF49f7hlRrQXXBYfWaDltTi4AgFpy1ZAQRQizdIGRTvFmRzAWZpeUPcLah6SFFdV5KH8CvivVwMfkVvyO6QwokWFv4lSvv1+SAIgZaSgl/iLdbimgjRATBDuCNMa0rdRhpSVEhKii1iIUTfGDj4kWYTcDPzxn6iWWd0n7FshBMS9mqG54NXUGFHtTd5bUmupDns8hRpFE0NtNGEYJzEc8eVQKSdGEVjvCAP19T9GPreKVbewXyWsq1PCBYp69kAVV1PAtM1QUur2omVt6gNAw5TpHvIbpWdhM6oUgtCFZEDO4SezhXwbJQB3MJ+HkR+tf9FAS2p0vWFCQ1ZQUbMqwm/zH0avyrbfCCNpYdWFBl88o3IVKltWCMQqI6OtmZ4Te1nFiZYReJ/Wtg9u+E3C6fomy1fmYzl6haCeRqBt2iZhemcyQpnpUB5vUEdblRhgiFGCQwWxRRY6epnhQtF4YxIsWUNSp7odjXxkwCfvFDwG36NWs0BfU6wPZTTdzVbevNfEmyDeUL7pu+ghyME58x1oMSvhp7KWU0qineD2WQrRCMJmPEANtmIxSceQ7If6x6s/5iWaegWgWURXRVfUNQi0u3MkqzoAqYUIwkLDrXtqQFYF1Mh06hREeXMtmBcrtPwPEQnJMMn/REAn48cJuof7ETdLMJm6RoEWqJ36pq3zZj0KCzi+QUSpKJKF2jWSfhU80BPUBzr92gSS+OkA9E2gUr9CGlctzTx4vdCxaQlus3kSrIdcyXSWmHuaqQO6dFiFdMjwhvlqtB/2a9ed/Nb2IwOvUAUPcrKI5z21pthZhiWWsBHl36hF+dgDjT94krknVtkUzJyEW7Ym9kdyUtqr5Zib/StwD2vWiNilcDhXAoeiJWsUKLQ21eZbOIPqmgXTO4CAZen4AAQ9PXrtSlLVJQSlgfRAVJn3DnnCyMn8AJiF0pp4QI1TUOEsH6LzEH5xgxFWlbDjQz7c2mbK7WPLvL4MIpCAiaZSxJILyqCpb1p4JSAFMCCnOILM7XlJFxzuQEznBJVqgG/k89msY+2iYSsH/5G+V6yH6V2BMP1xuGXFK/OBsB+wwEwRddA9OGXGTiaiSr83IgdwVUl0osseBlHLxKCjJ5S5C3MwT5+MSwqj1pIwe4iibOQUAwrGNJajaoy19KZA9IBTakymgeZGlXZOVRdzJ6F99m3J/A7hMJMJvwNlotuWRkR6U7w/DLALoYED0NASf4QCBx1ugBqSD/yn7RStewhdvproNy5VFQjuha8TRrJ0Id94e6Z9lj7lUTVrIJPIouxqTPQsARJ4jSOKU8maGkpVZnM0v0ytBwARs0c5ylYMkTi/dW2jb4qcPK/hpFdaDL2EYu/WKYAfHzEJBmPz3RLE4Fsd277buDwXG4dFoU4CFUztOlr7MSOR/lh2jphO0twWZzmYVo2bcKj/1dnImAA4HYEsQG7bT2D+Whk25VyiJnkjJGIRpk14cVEc1ZDnQdLdg+UE6hVOR6oOMqpeblKTrv2PCOcxGQ5U9NQWBWK+X9C5plkHq8xquFFpWpyItHKCncRPkHz4+LQHG7Xtii1JRahR7wa7VhN9Z9Z5HgzAT0REljQ9ZRQVyZYe+8vKgHvrRL4GV0YIxNZcNC21HhX6Uwwr0eAES8rPOo3QN3xqZBgF0/6U094AgFucFlhWxS/iZcpPTEfBHUTU+xfpUgwRMtKhXgywyY1CR/dUG0qK6tx2bGf1zBLk/T4HXgbZyPgKMM1K5YQ0G1ca2y7M1Xt17vdCGBbT/Y3mKpT7QWdYm61jSsQZ1VX4j1WUs7K5OluzghAW3P1hTVyrk0YB+4vFA8IE6FUkpUm7GahPRiWHpqJ2ZjqLeg7FqyZ+1qAp2lkAzbgQzBb0vAEQqakkpb6ezUpqUALxBqx9YwcWYSkewdga8fcNonlRBN6BfDor7ZhAaDkhMQScXYTGt7aR44JwFnMJBmeeWF5BoOgqVQ6K5i+L14pXIBaLoIlYxsJIaeFlQfpgxea4SkaAJrBGEm+6454AN4rpskUXl9WQiGXWLWA6WQrg3tUOws8LaZ7QfK21oNub7+Eeu7GtDkHyATGJ1qWbGitL0wvzgtAetCdnf72IJGbfFeqAHwhPtBl5qJF8clyK0ua+vUAy20OPko652MwzkDII8cc3lLjWYdOg9NzhLu4rwErNTpT7QhWVYElwF2TDW0Mo+Wam7ENJTOzGVOOTusDiexe8GF5pBKdRZ1QRfrU1gxnTZyy8ymKhVnJiBJpx2iuqQkJnJMpAF2T5jT+qqxfNxayERTLY1dmPp13j2Q29Gy104OSG0sq6mMLdtZlpMTcGQzghijl+MiQEZlktbCEWGtWmnK0lNMRUC1nPOKczVGpQ6lNYvco2XjqpabRRd8ewICmbZDUHGCjRLhLSHTBHMgO9CqyGYCvWREQIROgubYwTJJU9bWQq5RuzkxiinbW77zExD7g65g9TCkXEouetWECDhrciFqSK8LaxVpPASbgOSesncr0t3qmmMSU1abV6NaV82f8IB45n5JR3pc/R1WAHpmBWglZKh9eCKbD0BEG/DvQKoq1gBYtdFUWgxpoBRFT1TI3np4BwLi6OYTbdgmIVHuBNEC1l1BVcY2z6WbhFq3KJHVLK3kulch90bYYFn7XpxuL0umjHfNLiS9BQHH9yKWJCkXARuRz9gakHhUGSsGiQYh1qssqC3u9oL65qfVQR6FYlk8RmyzesO8Laia+8MTizbgOQSkKb5LVJckGrjOTDvglsaR96EF2NNYWBlYTROgUGaTQxxWd+BjPh2UYqtcw2UaOaA4H5XmhnchIJDUvysoTb3U5P1tlaQuDu9/9UiaO9ACq73eywuHcxv1gfrBYlNJtYfCQGKRh5B/Fywz6g9XOu7D7SLp9kQCIgOsFkw1uGJBwCwRa4coKQUTO2yRi1q6Rt4lN8LteBttuX3oNeVWIlVFBzGm2UKYnRRM9YA/159XIKAnDOv2oRoSK+ENadRBmVuN0rkxbUmWh8Qy1DnWW2oAU3QjHCmgsVGhmW0WzFNgCTmxsEgK09Qc8HJT6PZUAjo2uUwQRJqcy1PiJfUd9k1G8ORXAx6WstjYs9DlY7JqZ22s4OoyU1B6IqmwNCufHt1RAwBrN3kXfDgB+Xr1JDV3WZwblUQNcb6pBElSIx6/UwhYG1qOcskNAB25tXmFpid8ai8FRB5Q+lpy8TWMnxQZbVv8bgQcYCANaVIHk1R7IH1Rt0dzQLuViOOKhFBu5YZOdCI44DtBsjaPSxJQCLMQnCeE8HYEHKIgSlqkCr0BNBaP5z+AL/qyeCxPUJeXzqC+6I/BAc9hT+uiqdW24XFXzHwUuxNw8k9BeMDM2JFVHB5oliqawVw9FJ4AO76lI8cqK/6m5aT3AtQrh48J5SdaCBamaRl/PjlewgPCcCaY5NWMJdBryFmYGQPJRrVQ1SIM2SwbIuRkUs6vgdU0hYV4Slf4hkfD6rqEgvckoMzz25J2zCUasE9rLR2WKKTFa6i/7sGY3jigt48G7AFlpRZ/AdQckFYreFMCjjtBQ5CWo7jkIAFxI1Cd3Th/WBhuibjU+UbfyAENAmpKSKVsM/93QnbTPwjLEG1RWYNPZMoCJHvFUXCu7tKkLB+ojfoQlNjKDKNzCbX8X51/1+hUSs0a3+RxLAO6IWxR7SItaVp+j47RZCiVW9o+v0Do0sm1yG3ffvRHRCxi7ZJ0M5pbEHXfDAMEuX3d///4uvVFV+nfA7opTEnt0pYakKigoOIBlr9kg7wafFJr1XQG8tsN2Ec1OCgwa9l74K0e8Ovn99cIP/kPAndwKAFnZIJ4rWpR8YLm6qZMydK0LJ2+qGYJHZrV3gd8PdhCUB1WR6xl7zxgXCVbuAnyX/B2/x24QUYdTMBhBpJ4u1Qg54V11q2FWGMakDER+2jKAaH/OvRbL7M2eFpJ3DY31Qdbd8Ef98ddRhoM6t8NiiWGJavh+ftE+GjY38MYAP3eqxDrq+r11B4QoEhvE7DkczJnliNuGX0kBN9/j/rnXCH4DmoNhygoZk1J9YJoleSWxL3meiNdaEYS2FFR979GZzWfy2PSbJ2Fta1HxcAm5OMC+o9St/ACBBzfjDwuZEuFphOvLPECzD9uZY/oeo0Ka8Mkumo73Pp8xO+gzMRykgdciZcgIDtJaIui8y9LhSjRTlGoR+x5wF7Gj7t2iFmtHZ1k2gEiIJsr1A0y9Gi6CJn4GwRcEYelpCinOvkaA5burXp1nj0x8FK1110jZ62TESOHWgzocK9cYs2zaM3Yf4WAI06wmtJWIUoAOwulN3O5O2ANN2jiSk1FUBxcfVWqsgdEQ4GWut6i+B7J/0DwyKOWY+J7wmuVRTgf3OkaFFqhVnW5UidqyhCnLiK6/TLHA/Z0AN5l8ME8mkMlYPWAmqKGld/qY5lduO1CZNsqmmtKXcQcjOjyHz6yRsaMU/H8uSKfF5i6mpZOafxjmbfzekAYzwQBuAmVxSAl1fBi5awVYi/Vo4ypt0P6eExMTC8xW9A3enqAZKXJhzzg5XoPwNdPr/wDr0XAUQq6tsPldSrhlu6Aa7WyPEXI5i5tsT8BsYOzk2fyfnmj5KjUZtLEbgJefj91+XOqhxE0GObUZfF2UNWA7Ev/wcvZW26HEOrvGQQsZy2iJ4DCM0A+rUjrFhby1Zhuglwf94F/Lt4GD7wcAYfDsCIpVdTlUenR2HCO7iieEX4rldpXE9Tpgz71Wp2MjHDkXjBc/+eAV2+DQf1PhLSQR7RTjtZHWQdzta2KKWhod3QM+DEMRQHerRSfKXSwZ1uFhYc84PV8T8PoYDZaIamU67teJiqWUS020bjjNQkAOcwmO9hDlUbW4HaR4lsIePkB+Pw848MIGqjJV0nKCrRmqzYXDkwjX8s5wrLt5rOkrXFGl2M2N0tC2V/WzDD0RDRczvVEdAPU4l7RnorET2LQ/q9NCC9g1dmepcooRzvcjoY8bSjcUrU9akDbhfyZW3ECmiVbkiAElcVwrC9t5ZDZAY2x1UxCzpHKYUvg/2k/UFhXygj+LgFXxuGOhrKKyeWtRnbM4+j4r0ZFJwQr53maPsBUBt0H/mUC2sRqCXoo2MvK6Wobb3ZDPhLpCDVHLwwE0v/VyJstk/8g+I9hFrzBMQxCy6aaJAhJqYGW9NZ5LKebwtDlPM4haDtJ1UJ8LkICNnvA71M+Ed1Cz6xUUjnNkgpoQWeNnVzoYVSFh37KyIxETlpCldIKRwnyDgfRFA2D2ZKTKNg+bWs0cR/bDGjmLcEcmfGQvZSTIqJkmCBvcRBN0LSYT1RYHS8/NCKZkw2bmmwgoTZ04sjIVjiHYMDtwaLrgkGC3E72wXQfmM1XiPJlE4KNRU6tO15cekWwVRs0vS4OtFpLAPUzceXAWfZiG3R0E9I6iP65flx/1uo/FNRoXsk2BZ3LTZxll14rbsQZUbSvRplwPg3NqSsQHyi6A8iPStsGnUiQ7+v37/976d8VwvAeSTMOa5Ktpe6S4dkopJEVOLHIf4oJuB3qHwsTCXK/Wywe1zoLAZEj8klmExsVoqS/F7HeMrV2pd6g15taX6OrpfG3DrRtBR3tHA/o+lDS/VkZsUs+DQHRMUs7awbfXkQ4kNahmp8j4xngil0Loo06ghKClx4s89TXti3dBPn++r7B7YeHWKHqZb6gchxs4dZIigq2sF4KmBWDfOJDKPob32NQfJzeHeGUcSJNxjtpF7ycQH/bu+DTE9C+3WFJJiEplr9LwRUbWz/EkJKrN+PWWkK7IN0qpD+DnxRDT0S3G+gEPOBnGrYAP8HXkaTLbJczQvSXv8WM0V2w4f+a50KEZLwSy2kWqRLAdigSQ+RYzl9utgc8ew6YkRo21uUG9iLjxyiCHgMaVvWVt7dy7ItTs7YVrMsu/37hJsgj+/v+4gd9Fb+74M/LSv2vA77cHsEhCrLFHuGk0vEaNDc5LG2TQ0AjobMGEPlE19p+gnx/XT8+Gvz77/zu54DsoPqEBFxJwXa5ufww5hSltp68X3YZDulJzVqtvQcwAnqSwJkE+b5+iE3yKQnojsIDTnALK2pW2ntMz2w61hv9wBupa5Aqj5H03TP1334g1QQ1uVNyrRPky2++nxB8u93XemUi7Z0FVe3bygUBDTCjr5GUFSMcoEueSnYm06wVaCtoesDUsoe6b27DeScE3QxxNRjU/4rwW9GUFOVU0skWuWvxw5BubkJKT2LIuMq2CRttx8zhAW0wk/skN8VhkxjliHxaFG4AUj4r1AZcd026Kertt2bCmBEEbIAYf52kqKAFI5x4Yg7Ya1vSgZYdFvEO3vqnurbDbUrT6LK8wwSwKpHr6dPE1p8HYjZB3lb7xFTziCChFjD3Vtw5f6prM/iyjUuKlWehvbHAKj/a1S5hbWRLBbvZZmtWszveu8PA7/5TXZuhrZxDtFnOCjSCmCs/48kFdbuK9EPNAlsj0ExgWcw03fv/VNdmaPbsi7ZVMDoYK2wuvRetDx1bxzn906ISpJumSugwxxC/4w/8VNd2UOs7Bf1OkK95Iz9baLP9INCnRB9K0wPyHHYpNn3gn/ipru2gi+IT9O9F+Ltks2PwGFDdRnRb5OOWVn+tJLBuY/LVBro4xDGMF3SNnJIjFCTVUzK9DqwED+WArCYXoq/9s6de1AHM8YAr8S4EnOAEO5mghzMOGXcTc4eMHLLgFBR3lhQziDasQsPA41iX60d5LtWN9yHghExQqGAr38P+bvHRDbqZIetyYNfmzSeE6ixjDTyQ+nP/kRr+wGkPb0TAgRMZU1BU8MXts8Mhs6nBY1gAjQ8uYUkyOZFW9I3q/5Ly7+Us0NtgUP8p0DK/T1KWD5JlnhfU2bWMCbRaWOqdOeD0g+jy/wDei4D7O8E+bUYoCA0n1uBf3TUY9ca21tLWwqgH/H67r2cbhN++lqRcIlJg8GUqGixGY3y8Ui86wwMOXKAZgzngt/j2oVn6TwO/iRUgpVcAABDpSURBVE1JY8FzwQCPVkpBq5NlRMpYO9taS1UTIx9K+vi4jp5DvyEBJztBVdIm0fCxinLGWF6Zm1zHJPT7IAZfG4hzwBXw29iUFBU2ZQhtRu9sWLBcYN5nGENl7rEhiMWJEHsfBFyFtpE9knKdHPTIFPGW2roaPXRnoEr9VgAfiqqLlAQBV6K1ZKbg2ntzDoztjZsfy7OHjw8Je1NVdGlNvQT5zD9VffIfrJ6Hlp19krKclLQo5CrqqjFaWGNM+fZN/tOcqaJMa+okyOfn7fZ53wBfgoAZDUM3JDsVTe48EjStRmeSVdiJ13SESggG/QxGust5HvD+jRu3y/L9LyN4YwLO3oygEmgwyI9BH1qryPDq/ic7LzxabYJYlyZGSpwE+RX7gtZXwzQavi2axjYlzRqlQFtTlUJesoK4Z2tI4cFkx2jvQpIcuJEnsvcjBLyO8+/NCbjniQwM8cq/bVk6tcUBR0q9cX1Z5iBlp94J+RX7GOff2xNwpRNsb0ZoAXrdygCHCKgSJpVjwCoq764ldRdC9IGa/+kYImBP6nYRIm9PwKceS48cQzfSMq5WCmQfSBWVgYEkoBywz3ozCfhz/fmLBFzJwPUnMsNnhFRVK/5WhkPJFSkBAZGQzg3VqAI6vAR0/EwDXG5S3V8g4OwzQXszgkpdJJQe80E+sHciUDNAvrFIJfpa8+QaHJhMkD9KwN0zwYRyr0KHdQDA8VVWotEo48pNW5NTk0QLQcBJsHg1ICkqeAF1SW4fKE56ChGFLB6NHBdgYuozW1yr/qERDUHAWaBr6RZtljOdHsbxBrwVoCM6pUVztI+37YnVrNFntkkEybmhJODJfidkC7YzsOsEH/+OENCsA+kB20P1W8DrAXcgxx/2gHd0l1MTHHKC9jlgg2yCaPU4uk/APAZc0J0WID/bRhBwKtwMHNgOCwKO0A5F4bp/WY6bZXDWhgn5KxHqYDpz6u2VCYKAk7GKgh0NpMBNQDUFzDpAqbBHWXnXPGBRRtC110yCaAeFf46AkzNBraBBOkE0n2RjsES+d8CC1NctcRvxRPR0OFZVE22X0wKbWKxzRjbECicBEzvEabs1LOqZPwQBd8F2Bva2ww7yVcbQBlI/qpDDrVRKkDoHLLj3IOCh2E5BWS7p1SdgJi+io6I/SyrDZSrBuwtOaK/TRBBwH2xnYCcOdxlYGKP3oBIQBL9obwOzBsp3E0HAvbCdgk0GeigIhVeKg1Mb8G0D62do0q4GQcD9sJ2BytoPUrBQSttASBVMYHkAYqny3wop7RwtgoA7YjsFZfkgAx9NsBe01CfuIle4s0w/tTMdQcA94V8/U1BUKJwxuEcagJbiiQbM/9VjQN/DLY8eqj5HkyDgvvA7EDcDvU4w86kSkOtm8oKdgE9hXJPFLXyNgoA7Q/JnWLJJwbYTfNQClsSKeVs2oPK/+aigmAKUJ119KWAQcH9sp6AkCSloEFD9p6oV5yZyOCgK02IxeHRmAyKcmwgC7o+Gk/FKynJSYNAPp3GVTkZ3RTcWKH95K33kkHe//almBAGfATcD7ZxRqKAFOgWxA6tv9c7KiSHvmg1H39/SDYubfkHAJ0G4GZekWaNI9gioaqXlS/AUOR7NAVUCJoWA4QFfCm4GuiioSUqmWZxUOqpNQNCLjUSdAiIg/XBxD0HAZ0EnwZgkL28zUHN+hEmqACi7XCUfNcYsz3Y6CAI+D/6FsSV5ORXsUfCRApamotq3fdDPpdF+B+17ukYJAj4RA77BzUB06jsahqW74rtgAKVA34WgseAhdRkYBHwqVjHQfyBDKbW4O85FgPzlQwYxUWdZKr+ppdawCa+f/JmQI/SfDjtQkFJIcIo/0swOZhr0KxyiRyudsbcYLREEfDa2M1ByhkpyNqK7Y5APVZajFUpT3j/kY0N+PNNPElWFGoKAT8eAh7AkZXkStCMUhOq3SDGjKescCgFbiZ89bo9oEPAAbGdgOxWktYR+y+tyuNwaSqZeMo9nGgM/6MuJnq7/nPAz0EzRmgxEzx3QSiwnGqhdZxWeseK2XsEg4DHwU9AUFCpoQWZPoQ/yZaUFDdyyb/ONY8w+0SDgUdhOQamBFaT6P/Z4gLM59yiGp5bftGWDgIdhLgMtJygiLaZlQ/e2mRUn21UcBDwQbgpWRtkf6lBLxD6DvhgZw9C0io/t7l2CgIdiyAmy2GZo4AWJAUhoBlE5bVaPV1C4qGMqQT6vH1f2k+pBwDb8q888V0MDK2DuD9+RM44Mp0xqeQlP9ICXzxvc2M95BQE78DsgTB9ThVJAOAbleI/3PZGA6BiQ6ZTqZxLkl3q362763xRuBtJTk5YGTCl88PeoE1GZkHDenPKgtdKC6QQJAo5jgIKGoKAm8mkJhe9FdvnUWhJkHDnsa44TNG+nlU4nyCf9ResgoAd+BvoPBel7dksDuzzG0jnTyV3SYqV0NkF+vvbV/6agpPBKtisQsaCEv1ooOp3FP9uZ7ugB81dDf15kxZ/5nZBN8DtBS1KWU2LVmyOEjBzbZwLWfoaXTifH7fLJi4J8XqxiYE8FL2AxF6Tv3UzAh/czJrHvLhi+vkVRENCP7RSU5YKBYr/CCbidgUPBfOo5oPKL6kHAAXBv5JO0KtSSlB1UuVeW6/JDp7Nm4ZTe4XdCbnvp/wNY5QSHb4zgx+VRUjhEnObQvM/jQ9wLfjX4GTgch8l2uFaRXcr0Y+gugoCvht2cYCUaboB3w3On4JMNAr4ckiSQQ7JTkV/nrQiVzIeFcyfgEg4CviC2O0EZh1WiVYmJ/INym8+DIOArYp4TRMyS4Rc/ErPH8D2yQcDXxCwnuGx6UQHgZ6UG4+XY6F2iQcAXxTon2I/DgL3g1EcQRMcu0SDgy2KWEyTHywn937ipN2XkLtkg4OvC7wQbFCy3f5EgLIfFgG4O4yYTRh6bkLfABAaKp1Pw0/j5BaLLFHc44AKDgC8NPwMdmxF+JgiYgPiIsEd2z7i9njQI+OKY4QQTD7Tk265Is/6h4MDF0JH6RRDw1TGBgVoF3o1kf1X+bfbSE3gIeT9eEgR8fUylICmp7gwfCXYcIHRcZOmuLbMgCHgGrGKgeiLDfWBlYHGAjgzQeTDZVvRAEPAUmOUE+a24WlL+cQ2lLQHOQH1HEPAkcFPQFsxFhIH4QAZ8pOleAiP394KAZ8EEJwjY1yHR5VjaHTbdI3UIBwHPg1UU5DXsjA4FZvfRnX+cDuEg4IkwzEB5Q0ymem6dI8ME3ouJIOCp4KZgJaCkIFBXtwsD3Q41CHgupEEKJhBHwrIp+pqiuYN0yAYBz4bBOKwTjr+a7QXdpzBBwPPB7wTzbQtLDtUUdzltG+yUDQKeEKNO0BJMnIGTwrDf/wUBzwk/Ax/ZnU4spqIQkEqO8zHFNyO8P8YoaEhWyiFBenSzyiVGCP4D8GeCtqCoSIl7r5FoyvT6ZIOAp8VgJuhmID62cT18JbsboOD03wlhX9EWBNwPfidoUlBqYDLrPCC739dC/E7ImbGHE6S8GQimotHTD6J/vx81fqbhqRikoHYq3VGw6lTmoBB8RxDwufAz0Lzf5m0/PCqP5GyC3C4Rgp8LfyaYZewwPG9MIHvRMZkgHx+XXfUHFGwPw0MHJ84BPfdecPm9B+EB43dCnoBRCsrMbi4Bnep2IUfkgIdgLBPUdqjzMsChw8P4tcz3gD8TrGd7EzcdooNDbsXFOeCRmJAJzh2KT3YqQS7xi+kHYsAJVh+4DwUHNtVxL/iNMBCGe8+qbhwHuL8gMAj4TliVCe7CQPxQTbODIOB7YfhAZg8nmPBBIOuA9xUEfDMM+cA0lK4Nj2HpBbFO9hUEfDuMUNCfM44NoEZ2eigoT3+CgG+Iob3II1pOZiC6FYw5pxxRBwHfEaNOcHImSDomqsMD/hUMxNbpYZhra37GKQj4psgs8DJwOb2b17VVyd4HAd8W9SzEkwnOeyQrkRSwhyDgG8Of3xXyzfCBQ8eLQcC3Rr3d4X1OcI4LTId8Ku4I/YEekvqFG4rcvEe0alLZRxDw7eH+zvBZO+GhcB4EfHuMnMiIF6s6jMexAgSjB33bXOHIw1hBwD+BNOAEf+Vhkw9MiX3DUQtBwD+BIQpu3Q7XA0gHgoB/BNkteVPBbVmgP4YHAf8MkvdAZsZ22N08CPhnUMjnYNfm05ggYEAC3Zrb7TPBqCMXgoB/CgsBd/kkEu0mCBjQkPwPKGztRNMvy4KAfw05DdyRgtadYK3PIOCfw+ix9PouRDHEA6kB2OdzILIDyT/tgDsI+DfhPhNcq17THB4wkLF7GFa3IE/IAX+YwiDgiyLfsN37SJD1KkomE+R2DQKeBc84kOljMkG+voOA50Em35EMnEuQzx+uMAj4wqj3ho+j4FSCfH8JhUHAl0Yqv+t2FANnEuR2vQUBz4ZdM0GHypm/E3L5lArjd0JeHtkBTidgl9WzyfHxwCcpm6g/sA9KKjhbLzhYPZ0gEYJPiPIx3tnfEthndRAwADtlguEBA37s4ANdfI57wYGM+YfSz9sFH6Y/MBEzPaBXTxAwUDAvB/RrCgIGEGb6vyBg4CgMZJNBwMAOCA8YOBSRAwYORuyCA6dAEDBwKIKAgUMRBAwciiBg4FAEAQOHIggYOBRBwMChCAIGDkUQMLAD/E/VBAED0zHyXGEQMDAdIx9vCgIGZmPgN5mCgIEdMPLNg0HAwHTUL7/sywYBAzsgU7CPIGBgF0QOGDgYsQsOnABBwMChCAIGDsVMglx+v5/yupv+wBtiJkG+bvvqD7whXpqAWxpHx+/fMcc1CBgdP7Uxw/X76+P6M1H/Ka0SHT+vMcPX5w2+GQP/nlWi4+c1RmrK7z18X3lFINDAHAJiyk3XGAi48LsHuf9eXCBwBLQcMBB4Gm6Xq9gFBwKBQCAQCAQCgUAgEAgEAoFAILAnPrccTd8bf69se7usvyf9s2HQW/rdNOFNlv6Pn5Xj1j6H4cft6+N6WdnWg8vn7X8fK+2ypfHPda1B78/zfP/+/+x+N014k6X/43ZdOW7tMXh/r/9HfPu8bNDQwe8q3lZeHj8bGl9u65/Kudw7/rk8vd9NE95k6f/4+j6CgJ+f9797373dYJYtjVcT4feTBQf0+8AWa61v+/mzdtza5zDc+FqbcIzgdtnC8Mc1sgarifAxo/lqrJ/wBkvfn6FbS0DlcxhDjTe5UAc+Pi4bWv+sf7rwnATcMuHVlr7d3djaELzlGbyP/y7weycGloet11yXufGa/DS3PSUBNyXkqz3g5e51t4z7e2Xs/7iH4N0fYV6dmdwu68PRKXPAbRNeP+jl0xkbOl8dvzc09mDLvm5rirptF7zeFW0w54YJb7P0HSvHvelzGL8X3Fr36etgyzngtu35aiLcfs8BVycmG84BN0x46zngMTng7X/D702j7uGy4Wx/CQyrmLDpI3/f6we9qd8tE95k6Ufv65pt+xzGfRf8jKOYQCAQCAQCgUAgEAgEAoFAIBAIBAKBQCAQCAQCgUAgEAgEAoFzo/WoZ/yoRWAvkF9i+Sp/mNBThxT4SyDcupY/DaFAYCI+0IuP+6fGf//U7976uX8eIggY2AuYgI839z/lu7fuH7z7/goCBvaCTsCv/N1bv58GCw8Y2A1lD0II+JG/d+E3IbwFAQN7QfeAhJaxCQnsB52AV1IdHjCwG4wcMH8RwO+LDV/oGwi0QQj4+7UZj83vN9w+v+6nMPcvVQkCBvYCIeD39br8+bl+fFzuO5E4BwwEAoFAIBAIBAKBQCAQCAQCgUDgaPwD3zcfxCxPhDgAAAAASUVORK5CYII=", "text/plain": [ "# \"Fitted\", :ylabel => \"Residuals\", :term => [\"png\"]], @datasets=Hamster::Vector[#, @options=Hamster::Hash[:pointtype => 6, :notitle => true, :with => \"points\"]>], @cmd=\"plot \">" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "require 'gnuplotrb'\n", "include GnuplotRB\n", "\n", "x, y = model_fit.fitted, model_fit.residuals\n", "fitted_vs_residuals = Plot.new([[x,y], with: 'points', pointtype: 6, notitle: true],\n", " xlabel: 'Fitted', ylabel: 'Residuals')\n", "fitted_vs_residuals.term('png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Histogram of the residuals\n", "\n", "We can further analyze the validity of the linear mixed model somewhat, by looking at a histogram and checking if the residuals appear to be approximately normally distributed.\n", "The distribution looks in general not too different from a bell-shaped normal curve." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAMAAAACDyzWAAABNVBMVEX///8AAACgoKD/AAAAwAAAgP/AAP8A7u7AQADIyABBaeH/wCAAgEDAgP8wYICLAABAgAD/gP9//9SlKir//wBA4NAAAAAaGhozMzNNTU1mZmZ/f3+ZmZmzs7PAwMDMzMzl5eX////wMjKQ7pCt2ObwVfDg///u3YL/tsGv7u7/1wAA/wAAZAAA/38iiyIui1cAAP8AAIsZGXAAAIAAAM2HzusA////AP8AztH/FJP/f1DwgID/RQD6gHLplnrw5oy9t2u4hgv19dyggCD/pQDugu6UANPdoN2QUEBVay+AFACAFBSAQBSAQICAYMCAYP+AgAD/gED/oED/oGD/oHD/wMD//4D//8DNt57w//Cgts3B/8HNwLB8/0Cg/yC+vr6fn58fHx+/v79fX1/f398/Pz/Jf+lXIGH2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQ4ElEQVR4nO2dDXbiTBIE0Tm4z5xDP83e/wgLje2xwNs72apEqBzxnlnre0kn7gobIePZ0wkAAAAAAAAAAAAAAAAAAAAA4PcxTsM0fh6cvx8A+JmnuX5UlnM5lcu46wOC38UyXm/G5X5QPSzTbg8Gfh9TOT04h4DwQoZvt5WyjLs8EPidPAo4DMteDwV+I//wE3AAeCZIwH84B+yq6nt8Kave/fG9suqZZbzenD+edm+fI+Bb3Clp1TOlXgcs9xV/vg6YcysQcI+qH5in4X4d+rbi8nmwuerttwIB96jq4oVVcBgQEHYFAWFXEBB2BQFhVxAQdgUBYVcQEHYFAWFXEBB2BQFhVxAQdgUBYVcQEHYFAWFXEBB2BQFhVxAQdgUBYVcQEHYFAWFXEBB2BQH9/BHZ+/G+FAT08+c/Egh4/Kr3AgEbIKAfBGyAgH4QsAEC+kHABgjoBwEbIKAfBGyAgH4QsAEC+kHABgjoBwEbIKAfBGyAgH4QsAEC+kHABgjoBwEbIKAfBGyAgH4QsAEC+kHABgjoBwEbIKAfBGyAgH4QsAEC+kHABgjoBwEbIKAfBGyAgH4QsAEC+kHABgjoBwEbIKAfBGyAgH4QsAEC+kHABgjoBwEbIKAfBGyAgH4QsAEC+kHABgjoBwEbIKAfBGyAgH4QsMFWK8ZpmMbPg/P1YL7+7zLcmIKrDgsCNthoxXwVbq7SXVnO5VQu4+l0KYaq44KADTZasYzXm3G5H9w+P5UJAR9AwAYbrZhuqpXvT7a3gwkBv4OADTZaMTwtcj5fBZwv384Mg6qOCwI2iBZwvFxvLteTwfnJQAREwGeCBTwvX5/OT6+Cb2yrOyYI+D8IMGJ9DliW8/flH+u2VR0XBGwQ8Sr468fe5eN6TH0NMl9iq44LAjbYaEWp1wHLfaFl/PivnAOuQMAGW62YP375cVtouFOuT8UTr4L/goAN+F2wHwRsgIB+ELABAvpBwAYI6AcBGyCgHwRsgIB+ELABAvpBwAYI6AcBGyCgHwRsgIB+ELABAvpBwAYI6AcBGyCgHwRsgIB+ELABAvpBwAYI6AcBGyCgHwRsgIB+ELABAvpBwAYI6AcBGyCgHwRsgIB+VAE19v7qNoKAflQBtfTeX91GENAPAjZAQD8I2AAB/SBgAwT0g4ANENAPAjZAQD8I2AAB/SBgAwT0g4ANENAPAjZAQD8I2AAB/SBgAwT0g4ANENAPAjZAQD8I2AAB/SBgAwT0g4ANENAPAjZAQD8I2AAB/SBgAwT0g4ANENAPAjZAQD8I2AAB/SBgAwT0g4ANENAPAjZAQD8I2AAB/SBgAwT0g4ANENAPAjZAQD8I2AAB/SBgAwT0g4ANENAPAjZAQD8I2AAB/SBgAwT0g4ANwqwYp2EaPw/O14PZVnU0ELBBlBXzVbj5U7rlXE7lMpqqDgcCNoiyYhmvN+NyP7h9fiqTqepwIGCDKCumcnpwDgE/QcAGUVYMT6udz6aqw4GADWwCjhdX1Rsg/r9pIeD/xiXgeXmO3Aiq2xnRESn9iwSMNGJ9DliWx+ffU66fgEajfpGAN0JfBX/92Ls8XQQMrHoDREekNAJ2Uep1wHJfcRmdVW+A6IiURsA+5s9ffgwfz+3DUExV+yM6IqURMEGVG9ERKY2ACarciI5IaQRMUOVGdERKI2CCKjeiI1IaARNUuREdkdIImKDKjeiIlEbABFVuREekNAImqHIjOiKlETBBlRvRESmNgAmq3IiOSGkETFDlRnRESiNggio3oiNSGgETVLkRHZHSCJigyo3oiJRGwARVbkRHpDQCJqhyIzoipREwQZUb0REpjYAJqtyIjkhpBExQ5UZ0REojYIIqN6IjUhoBE1S5ER2R0giYoMqN6IiURsAEVW5ER6Q0AiaociM6IqURMEGVG9ERKY2ACarciI5IaQRMUOVGdERKI2CCKjeiI1IaARNUuREdkdIImKDKjeiIlEbABFVuREekNAImqHIjOiKlETBBlRvRESmNgAmq3IiOSGkETFDlRnRESiNggio3oiNSGgETVLkRHZHSCJigyo3oiJRGwARVbkRHpDQCJqhyIzoipREwQZUb0REpjYAJqtyIjkhpBExQ5UZ0REojYIIqN6IjUhoBE1S5ER2R0giYoMqN6IiURsAEVW5ER6Q0AiaociM6IqURMEGVG9ERKY2ACarciI5IaQRMUOVGdERKI2CCKjeiI1IaARNUuREdkdIImKDKjeiIlEbABFVuREekNAImqHIjOiKlETBBlRvRESmNgAmq3IiOSGkETFDlRnRESiNggio3oiNSGgElxmmYxs+DstyXW4YbU3DVGyE6IqURUGGe5vpRGafxvtylGKreCdERKY2ACst4vRmXj4NyQsDNRiGgwnRTrfx9sh3+/tfoqndCdERKI6B+9+HheJov384Mg6reCdERKY2A+t0fBbycy/XMcIyteidER6Q0Aup3fxSwMj+9Cr6xre5dEB2R0r9IwAAjfj4HfP78p+MDIzoipX+RgDciXgWfl4fl6muQ+RJb9U6IjkhpBFQo9TpgOa2fjDkH3GAUAkrM03C/Dj18PKXfntTLMvEquNcoBExQ5UZ0REojYIIqN6IjUhoBE1S5ER2R0giYoMqN6IiURsAEVW5ER6Q0AiaociM6IqURMEGVG9ERKY2ACarciI5IaQRMUOVGdERKI2CCKjeiI1IaARNUuREdkdIImKDKjeiIlEbABFVuREekNAImqHIjOiKlETBBlRvRESmNgAmq3IiOSGkETFDlRnRESiNggio3oiNSGgETVLkRHZHSCJigyo3oiJRGwARVbkRHpDQCJqhyIzoipREwQZUb0REpjYAJqtyIjkhpBExQ5UZ0REojYIIqN6IjUhoBE1S5ER2R0giYoMqN6IiURsAEVW5ER6Q0AiaociM6IqURMEGVG9ERKY2ACarciI5IaQRMUOVGdERKI2CCKjeiI1IaARNUuREdkdIImKDKjeiIlEbABFVuREekNAImqHIjOiKlETBBlRvRESmNgAmq3IiOSGkETFDlRnRESiNggio3oiNSGgETVLkRHZHSCJigyo3oiJRGwARVbkRHpDQCJqhyIzoipREwQZUb0REpjYAJqtyIjkhpBExQ5UZ0REojYIIqN6IjUhoBE1S5ER2R0giYoMqN6IiURsAEVW5ER6Q0AiaociM6IqURMEGVG9ERKY2ACarciI5IaQRMUOVGdERKI2CCKjeiI1IaARNUuREdkdII2Mc4DdP4eVCWH9ZFQEccAe/M01w/KuM0ImCvUQjYxTJeb8bl46D8tC4COuIIeGcq15sytdZFQEccAb+vMzwcW6reANERKY2AG9ZBwO1GIeCGddoC3giq2xnRESn9iwSMNIJzwCijfpGAN0JfBZ+X1rrvLOAfEaNRCNhFqdcBy+n5yTi8yoHXEevie2/dRsKsmKfhfh16+Hhuf3p2R0DL4ntv3Ub4XXDF64h18b23biMIWPE6Yl18763bCAJWvI5YF9976zaCgBWvI9bFNfbe6CcQsOJ15H0WR8A3RRsjAsaBgBVtjAgYBwJWtDEiYBwIWNHGiIBxIGBFGyMCxoGAFW2MCBgHAla0MSJgHAhY0caIgHEgYEUbIwLGgYAVbYwIGAcCVrQxImAcCFjRxoiAcSBgRRsjAsaBgBVtjAgYBwJWtDEiYBwIWNHGiIBxIGBFGyMCxoGAFW2MCBgHAla0MSJgHAhY0caIgHEgYEUbIwLGgYAVbYwIGAcCVrQxImAcCFjRxoiAcSBgRRsjAsaBgBVtjAgYBwJWtDEiYBwIWNHGiIBxIGBFGyMCxoGAFW2MCBgHAla0MSJgHAhY0caIgHEgYEUbIwLGgYAVbYwIGAcCVrQxImAcCFjRxoiAcSBgRRsjAsaBgBVtjAgYBwJWtDEiYBwIWNHGiIBxIGBFGyMCxoGAFW2MCBgHAla0MSJgHAhY0caIgHEgYEUbIwLGgYAVbYwIGAcCVrQxImAcCFjRxoiAcSBgRRsjAsaBgBVtjAgYBwJWtDEiYBwIWNHGiIBxIGBFGyMCxoGAFW2MCBjHVivGaZjGx4NluDEFVznRxoiAcWy0Yp7m+rE+uBRDlRVtjAgYx0YrlvF6My4PBwi4IY2ACtNNtTI9HEwI2J9GQP3uw8PBNF++nRkGVVnRxoiAcVgEvJzL9WRwjK2yoo0RAeOwCFiZn14F39hWZ0MbIwLGEGDEz+eAH8s/1m2rsqKNEQHjiHgVfF4eDuprkPkSW2VFGyMCxrHRilIv/ZX7Ql8HnANuSCOgxDwN9+vQw7eDsky8Cu5NI2CCKhltjAgYBwJWtDEiYBxpBfyjIY0RAePIK6BxjAgYBwJ2jBEB40DAjjEiYBwI2DFGBIwDATvGiIBxIGDHGBEwDgTsGCMCxoGAHWNEwDgQsGOMCBgHAnaMEQHjQMCOMSJgHAjYMUYEjAMBO8aIgHEgYMcYETAOBOwYIwLGgYAdY0TAOBCwY4wIGAcCdowRAeNAwI4xImAcCNgxRgSMAwE7xoiAcSBgxxgRMA4E7BgjAsaBgB1jRMA4ELBjjAgYBwJ2jBEB40DAjjEiYBwI2DFGBIwDATvGeGABRfxzQsCOMR5YQCn9ih+YCNg1mKMurqURcAPewRx1cS2NgBvwDuaoi2tpBFzh/DdP1cEcdXEtjYArnFutpY+7uJZGwBXOrdbSx11cSyPgCudWa+njLq6lEXCFc6u19HEX19IIuMK51Vr6uItraQRc4dxqLX3cxbU0Aq5wbrWWPu7iWhoBVzi3Wksfd3EtjYArnFutpY+7uJZGwBXOrdbSx11cSyPgCudWa+njLq6lEXCFc6u19HEX19IIuMK51Vr6uItraQRc4dxqLX3cxbU0Aq5wbrWWPu7iWhoBVzi3Wksfd3EtjYArnFutpY+7uJZGwBXOrdbSx11cSyPgCudWa+njLq6lEXCFc6u19HEX19IIuMK51Vr6uItraQRc4dxqLX3cxbX0gQUcp2Eag6ucW62lj7u4lj6ugPM014/QKudWa+njLq6ljyvgMl5vxiWg6tudnFutpY+7uJZW/zGtjQMOZCrXmzIFVCHgjml18Y0DDmT4ae2nKvH7S9wNY/q4i2vp9AJad8OYPu7iWjqbgADPWAT88RwQ4FUs4/XmvPyfFICJUq8Dlr0fBvxa5ml4vA4NAAAAAAAAAAAAAOBkqb957vj18Kj/uvr8wx8E/NOd9GvnZdEf3k9/r2Ap6vyiurbvhjyqbik6uHT+Yq5M8sYv53Iql/EFd7q6pH9//Pj3Co6il+7EqWdUvVL00Nt1mfUJ327U9+KMPXc6LUV/Q9tyq3r8ewVHUecX1bV9N/RRvVLAzrcmnMfOtyx2vRms507yw+t9r1rvO+detRMdo3rl+1Wm+dJxZjFfOje+LHLXlfNZv4/88H58t66j6IOeL6pj+3pG1SdFH5frmcWslpXbt0jPxg9Dz5sRx0tPVd8dXiZgzxfVsX1do+qSQufv+67nf//BXu+03L55la/qq0r5Fv68k/YO2s97vbmAfW8L1n8CyqP6QpBiO+Ij/PizgY7nEP0kpiw9NW9+Dtj7RemPb8OoXvPPw9TTzfklz3Fj7VN38NL5Btq+V8EdP5h65tTzRY23m74/6FEfYr8UOv1P9y+6Dijmv5AfXu/fK/RcBxz1+/RfB3zbc8BKWXqvr+sbv3Rc//94ChG96PrTwp6/V+j7G8a+L6pn+z76xHy/FAAAAAAAAAAAAAAAAAAAAAAAACDxX66uA7frdpvBAAAAAElFTkSuQmCC", "text/plain": [ "# \"fill solid 0.5\", :term => [\"png\"]], @datasets=Hamster::Vector[#, @options=Hamster::Hash[:notitle => true, :with => \"boxes\"]>], @cmd=\"plot \">" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bin_width = (y.max - y.min)/10.0\n", "bins = (y.min..y.max).step(bin_width).to_a\n", "rel_freq = Array.new(bins.length-1){0.0}\n", "y.each do |r|\n", " 0.upto(bins.length-2) do |i|\n", " if r >= bins[i] && r < bins[i+1] then\n", " rel_freq[i] += 1.0/y.length\n", " end\n", " end\n", "end\n", "bins_center = bins[0...-1].map { |b| b + bin_width/2.0 }\n", " \n", "residuals_hist = Plot.new([[bins_center, rel_freq], with: 'boxes', notitle: true],\n", " style: 'fill solid 0.5')\n", "residuals_hist.term('png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Q-Q plot of the residuals\n", "\n", "The Q-Q scatter plot deviates a little at both ends from the diagonal. However, it isn't horrible considering that we are modeling count data with a normal distribution." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAMAAAACDyzWAAABNVBMVEX///8AAACgoKD/AAAAwAAAgP/AAP8A7u7AQADIyABBaeH/wCAAgEDAgP8wYICLAABAgAD/gP9//9SlKir//wBA4NAAAAAaGhozMzNNTU1mZmZ/f3+ZmZmzs7PAwMDMzMzl5eX////wMjKQ7pCt2ObwVfDg///u3YL/tsGv7u7/1wAA/wAAZAAA/38iiyIui1cAAP8AAIsZGXAAAIAAAM2HzusA////AP8AztH/FJP/f1DwgID/RQD6gHLplnrw5oy9t2u4hgv19dyggCD/pQDugu6UANPdoN2QUEBVay+AFACAFBSAQBSAQICAYMCAYP+AgAD/gED/oED/oGD/oHD/wMD//4D//8DNt57w//Cgts3B/8HNwLB8/0Cg/yC+vr5fX18fHx8/Pz/f39+/v7+fn58AnnNucONpAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAV6UlEQVR4nO2dCYKjKhQA9Rzex3OIS9//CC1uETUJEBCQqv87M230SZKaxxLEogAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAKeUn5/uxH6HLzu/CXx5lEkoSI+hLct2+LRhQRGhPT/d75+yswYB86NpR3H6tnm/YUURQbx/Wpx21gYBs2NYUtmW8k4byr4thfxlEmEQ0y9lWTbrAdOWftwymzI/Vfbj5u51wIzoZT29bupkqu3WwG3Z9uXyi3zomnG3Zv1l2RWexurZ6t15w+hF0U/Sya2iHxt7/S4tnbfMCo26VGOIXvTTz/JU0+02yVPJ08gD+tHW0fSdgO3QjRI26y+7EsGDEEta6cS7DaU0sl88kdXz9svEecus0Li5Kxefh2Z9Su68bSp3B8xh9gJupSj3G+BhlMe/nDdIIzcRpl/K/W7nLfuH7ejdU+umphm6deu0R3cSsFi2rLvC09gnvKkVd86A62N5/OX89GvLJuCM8tS2qRJbG1ANMz1UQky7Lb/QBnwkGm3AfQYU+18mzlsUAdXe8rR1v6lbK+/pGSUDNlXX7c2edoWnMSw6vHrBxw3lvpGn3wZcHtq1//F66rRpCzwsdXa/KdmXx8DwML6OA8694P7VC566xFvHdtvyyonyqVVAeWRXbdlUPmyb9r3gYaxgpzGatpJjLjLMUHRDOyZYesHPZhDHb0LUDZfjgKMs4rX76+lifWpLgzJas7beyu0E03jM1TjgaJ/oZWx52qarynmUhnHAB9OX3YcNVHzgm+qYXPYbEBC8U5Xi7QYEBAAAAAAAAAAAAAAAAAAAAAAAgEfSTFcViu87AviAaxMgKAgIQREICCER63WKACFoq26+FBsgGP2hF1wCfMS1geXHXy2jeD2KkyV4siumPshxqaYUXgwnC3WyP6cCXrYBH/rOcTIHh/392Z/simn17GMf5JHvHCdzcdjfLyfThgUs4JK/2T8EhCD8rX9BQLifv80/BITb2emHgHA7f8pvCAj3ovqHgHArfwf/EBDu5KgfAsKdnP1DQLiNU/UrQUC4iSv9EBBu4jL9FQgI9/BGPwSEO3iX/goEBP980A8BwTuf9ENA8MzH9FcgIPjli34ICD75lv4KBASPfNcPAcEbGumvQEDwhJ5+CAh+0NQPAcEHuumvQEDwgL5+zgUZjvEQMDsM0l/hWpBOIGDuGOnnWpC2R8DMMfTPrSDVcIqHgFlhVv0WtVtB5NqACJgzhvpJXC7PJpfJR8B8MUx/Mv8VLgVsquJCQB/LAEOEmOo3JsBRDZcCLotOV+pGd/EhYoz1K9xnwAmq4Cyx0M99G3ACAXPEyr/CeS9YgoD5YVP9bvBdMPzIF/3q5Wd+2P5YQUD4jTf+1fP/KsuG/X4ICL9wWf3WhSqd3LQkwi0friAg2HOt36be9rgKWNebjCsICNZc176zcrsUSAYEH7zp/Cq579AIpA0IzvjU+CtU6+gFg2uu0t8p251suwABwYLL2nef/7QjISCYc5X+1KEXbRAQTLmofk8DztogIBiiMfRsAAKCEV97H4bxEBAMOOmnqGeuHwKCAUf96jPGMREQdPGgHwKCNp9qX0v7CgQETb5Xv3ZxERB0+Jb+LPVDQNDic/r7JTICwlfU6teZexMICN94r9/v/iEgfEFJf671Q0D4wnv9nPiHgPCJffpz1vFVcClIJUoxeIwPN+NfP6eCNFVXdO3BQARMlw/6uTuJQ0F6+dAJb/HhVnbpz5t9hXtBEPAh3KOfc0G6hir4CdyU/grXgpRl4zU+3MNt+t2QAVkjOjle6c+zfT7koA2YOpt+/nq+exwKMuU+BEyce/VjHBBUFv/u0s+tII0oRe8xPvjm79o/n6fku2DYuF8/BIQXs3+36oeAsDJXvzfrh4CwcKHfLedFQCiW9BdAPwQEyVm/u/xDQAiY/goEhKD6IWD2nDsf954fAbPmYuzl5hIgYM4E1w8Bc+bU+gtRCATMllP6C1IKBMyUc+c3TDkQMEsCj73sQMAciUY/BMwRmf4i0Q8BMySi9FcgYHZElf4KBMyNyPRDwMw4+he6PAiYFcfqN3R5JAiYD/GlvwIBMyLC9FcgYDYcqt/QxdlAwCyIsfU343qRcpbmiJEoW38zLE70fOJNfwXLs2VADNNO38MClQ8n6vRXuBekqvzGBzMibv3NOBZkaP3GBzP+ItfPtSBVc4rPIuXh+Puro/bPtRxdU522IV84DvpF59+ES0Ha4yCg4/hgQuzpb8HlOODxTpmO44MBiehnIEjXyp+y7d6HmlF3QMAQ/MXf+VjRFkR+xTHWsdWxm+sqPrgjHf0MBBl37GR2MzQKAe/nLx39zAQchMkBhvHBEWrrL3RpvmJSBctRltNIs6v44Iak0l9h1AmZbsV6NdTiJD644C81/5iQ+iTS0w8BH0SK+pkI0jdi3PnDOOCP8eE3ktTPQJBBDF15Nd3AUXz4jbT6vi+0BZFXe0xjgZ7iwy8kNvayw2QccPvxER/sqZPVzzwD9uL7rlbxwZqU/TNtA/biasqLi/hgSa1Uv6FLY4xBL7gtS2E6Do2AnqnrpPVjHDBxlPQXujBWIGDCKOkvdGEsQcBkqdNPf4WuIOUOH/HBmGfoRwZMlPoJte8EAqbIc/xDwASpleo3dGl+RLMNuGsG+ogP+tRK+gtdmp8hAybGo9JfYSBIs/zJd8EheZp+hrNhCiYjBOV5+mkL0oqlBXhaBHpH15yjIaAzHtb4WzDNgB8YxICA/lD8C10Yd7hcnKi7iIaAbqiV6jd0aRziVhAE9MRD019hdmH693FABPRCvfcvdGEco780R6VxQSYCeqDeV7+hC+Mch52Qy51YI/pH6genPyM5hM4V6WRA1zw7/RUGgvQ6l4MgoFuenP4W9DOgzmQEBHRJ/fj0VzgV5FJQBLQlC/3MBWFpjnuoldo3dGk8oi9IM1XC4nwvGkfxYUcm6a8wmY41fdM2sDzbDWST/gqTTsi8Qv7QeIoPK7WS/kKXxjdGA9HifD9gZ/FhplbSX+jS+McoAwruE+KZWvEvdGnuQL8NOEx3o+Y2DT6pleo3dGluwWA2zPjQfLpX3G/xIcP0V3BVXDzUSvoLXZrbQMBIyFM/BIyEXPUzGYZZL4vzFD9rcmz8LRgKojUp64f4WaJ0PkIX5m6MJyOQAR2jdD5CF+Z+jAVhINopddbprzAWpKsYiHaJ4l/owgTBtBPCQLRD6syrXwnDMMFAPwkCBqKm9p0wHgc0vM4XAS+plfQXujQhMbkssyu64dPybD/Fz4ma9LehvzTHbF5PL/hnSH87jJfmYBzwV5T0F7owwdHPgPP4CxnwR5TqN3RhIkC/DShkG5D7Bf+IUv2GLkwM6AsyyAUC20/+DaI8+YmAe0h/JxwKMubI6cdX/PQh/Z1xuUa0zH7H64YR8AX6XeBQkGkFweN0LQRcQb9LXK6OdRUQAWcYe3kDAt4CnY93eBeQNaLpfLzDsRy0Aa+pqX4/oCeI1lSYqRdcNTbxnwzp7yMOZ8N00zjgYcZ07gKi3xdczobpxflmmpkLyLyXbzAbxic0/r5iOhumYzaMNoy9aKAtyNz66z/ORvgl/vOg9aeDviC9vGu6qX/5Ckj604Or4vyAf5ogoBeofnUxqIIbUW59EQ/xHwTpTx+DTsggb9N1/KbDWfwHQfozQP82Df20M/eK+wbpzwijgej1x0f8p0D6M8M0A/YsUPkRvnozxLANyGWZH+GrN2NMBqLLUpguEZ2VgDX+mcM4oDNo/dmgf6suz/GTh/RnhX4v2PhrYLP4qUPvww79DFi1ZWPuYCYCop8tRoIMjWhYoPIC9LPG+E5JDESfeTX+8M8UwwxYispj/DSh8/ELRm3AtjLuCz9ewH3rDwHNMegFm9tnEj9Raqrf3zC+Ks5X/CSp8e9X9Ccj2I1EP1nAmur3d/RXRjAdgDGMnx6kPxfoZ0CN2yR1zenZ5wrI2LMTXAoyiCEfAUl/bnApSNOdwz1UQNKfKxxfFZeHgOjnDsdXxWUhILWvQxxfFZeDgIy9uMTRVXFr7/gs4OPWiCb9OcNIDr2r4h6fAUl/jnF8VdzTBST9ucbxVXEPFxD/nONSkKuvSh4lINWve7gsU59X+sM/Z+hPSG3lT5nv8mykPy/oL1I+TLdqqHJdpPylH/65xGQcsCu7bFfHeumHf04xEXAQJgcYxo+brfGHf44xqYKbapQwyyoY/7xh0Akpm2K7Y5f7+DFD9esPhmG+Q/rzCAJ+Bf98oi/IIL+KM74yKX0Bt+oX/3xgdq+4Lrslel/64Z8XTK8L7vJapHyrffHPE8YrI+Q1DrilP/zzhP444Pc7pv8UP0Ze6Q//fKG/MoIYurENmNM44Cv94Z839AQpd/iIHyOkv1tgHPCaGv/uAQEvqXfVLwL6BAGv2KU/9POL/kB0Kyzu1JWmgKS/+9AVpGllJ3hoG0/xYwL9bkRTkG0aoOn9ktITEP1uRVOQre59/ED0Sz/8uwNNQcqLv7mMHw34dzMIuOcP/+7GaRVcjR3lQz85KQH/dmv/hS5LLrjshDTV2FE+7JCQgH87/0KXJR/0h2F6OQwjmg/7TOodJgwmI+B0l3Py3+0YDkR/H4RJU8BJP/wLgHNBKvV2mmkIqOqHfzfiWpDjhetJCKjUvgh4K44EWecJnhbRT2CN6D/8C4RzObrmdDvryOUrTtUv/t2LU0EuJstELyD+hcXprbou+siRC/h39C90gbLDwxrRyhqqcQs46Uf3NyQ5z4gm/UVAxgIe0x/+hSBbAf/wLwoyFXDRD/+Ck6eAi374F54sBTz7F7Q4WZOhgH/4FxH5CUj1GxXZCXjufeBfSDIT8Dz4gn9hyUvAi/SHf2HJScCr9Id/gclIwIveB/oFJxsBLwZf8C8CMhFw04/Rl8jIQ8BNPwafYyMHAa/SX8jywI4MBLxKfwgYC48X8JX+SIAx8nQBd/rhX4w8W8B9+mPplyh5soCqfiy9FiUPFlDVD//i5LECHtIfAkbKUwV8px/+RcYzBXyb/hAwNpwKIhcpP6wPE0TAo374Fy9OFyeKZJHy9+kP/6LDpSDT4mzB14j+UP0iYHw4FyS0gKfql/o3alwL0jVhq+APzT/8ixHHgpRl4zX+F47Vb031GzuOFyk/Z8AbFyl/880v+sWKFznCtQGvrnnDv9hxKUjIW3XV+xtdImA6PGEcUAqGfoniVJDm9tu1LoJ9Tn9+iwA/kfJ3wZth6Jcu6Qr40u+Df97ODo5IVMCXYu/183NmcEuSAu4sO+vn44TgjbQF3KU/H+cB/6Qo4Cn9eTgH3ESCAirpz314uJX0BNylP+ex4XZSE3BX/TqODEFITMBX9es2LoQiKQFJf88jJQH3vQ+HYSEk6QhI5/eRJCMg+j2TxAT8Y+D5YaQi4M4/RxEhCtIQkOr3sSQkIPo9kRQErPHvuSQgINXvk4lcQHXilZsiQUzELODrmzd6v48lXgHrWk1/CPhIYhWwPqQ/9HsocQqIftngWsDhENAmfr3zD/0ejmMBO/GzgPXRP0dFgyhxLGDbuxDw9dWbo2JBtLgVsBqOAc3i10r167BgECtOBezbU0CT+LVS++JfFrgUsBPdDwKiX5a4XCO6qc4B9ZcBRr/scL1GdDlTKdt0DqyXrsc278VhoSBunA9Em1fB9QHSX07EIOBu6KWg+s2M8AIu7s0PF3fagkcT/rvg+tUCxL/8CC0ggy+ZE1zAYsmB43/olyGBBaxXAS/u8ws5EEUGlH9BvzwJLuDS+cW/TAkt4JT/0C9fwgtI7Zs14QUk/WVNcAHRL28CC0j6y52wAqJf9oQUkPQHAQVEPyhCCTjNe9m+BIGMuVnASbia+6zCyg0CvhRbhJPzXuo1/yFg3ngXcJ/lZuGmaX/rNNSCFJg3/gUstiy3zL3626ZgkQHhHgFXA+XP3+EyOM8FgLi5MwMW26z7en1Av9y5tQ3IVx9w5NZeMP7BkTvHAfnuA07cJyD6wQUuBWymtYnEdXz0gytcCth2b+OT/uCaewQ01s+uWJYvhpOld7JLxBsBLdLfc985TubkZJeIvi3FcIpvU/s+953jZE5OdklbdUV/MDCFF8PJ0jvZKcy21m8vjk8AfMCNgHvlnEcE0GLqg8h7hQCE4KoNCHAbXSNOvWAAAAAAAAAAAAAAAI9cTdLXZDD/TrmyGweXh/XGR3WNzZfeg1UR7c5l+cIs38XC5iP7wQ89ruZI69EJ81dTdUXXGr93docNwuJfSNGPQvTGUtidy/KFWb6LhdVHZu+H9xO0vfGrmT7Xzvif02B1WNPZTPtp5MmG5pZzWb4wy3exsPrIvAt4NUlfi2qwnNZl9dbZHWZRwOn9uOlcM3bvh81RNh+ZtR/aJ7iYpK+DnNNl9Z53jV3zparMj7EoYPnbgRbYvDCrd9HqI7P1QxvLCVqd/Jdh856XZWNx1Fhb2UxiTEFAuxdm8S7afWQ+J/CV7ybpaxzWyH+2Jq9mO5nZv931sKoxOGg7KgEBzV7YhkUGNP7IXhj5YYVxuZarBaxqD/NX0zU2J0qgDWj7wmyK+MNH5vMijl8m6RsXy67XJ3tvxodM2PaCbdKS1Wdk9cJs30WJcSm9X8TxSx1/2zigbRPEQopuGge06PpZjQNavTD7ccDI2oATv0zSN3/PG6uR/6XyMLTC8prC3qaIlueye2GW7+J8RtMDuIgDAAAAAAAAAAAAAAAAAAAAAAAAAPKmfPuL2aF2Ieb9uSfGkyiXicCan+q2W6t/zLzzl3gGURDwSZTLShGmAgr9Y+adv8QziIKAT6Jcrs+Sn+q6jJXoRTk9iH58kIZO1zE0xevDL+WF3KXcZfJ3WwBLCbH+Nu08Hzq8iVf0bTmts1Kuhdl2EUsZ5lMWr0DTibq2LL2vtgL+KJeFLCZd+ukStqJsujk1VkI+jIa2Qzca0eyEmUxYnn0tgKWG2DZvVvXjEcNVPHng6OBu19cu61nK3c8WWl7LZrUUAsRBuVxKu/45LXcybRgfuvlh2VVeOKsKuD7brgtgqSG2zS+r9pee7ePNB+4FfO2ynmX/o4aGdCllLTYvcTJdzNgpEqg6lEcBt7+six8cQmyb1x3Vy4N38aY/u5OAxX6L8rOGbpqBCjhl5opxJ8O1gJUQ01W6bwRcL+I9hFA2KyMox3gX51Z2OQm4XTZcCdqAKTN/1pX8U5wS1vbQVF33IQNuvVw1hNjvrGTAU7xz9lV3OQq4hZaHcbPShJkNmDqgrwbc9sT6MP3ZvxVwa9qpIbbNpx1P8aYnhqUd0JenXU5tQGUdAxqC6TJ/dt08aNLPi5acBRy3dkMrdsvmqp1b2fmVPVU1xLZZ9NuOcy/4FG8QcsEU+a+gkiMrxXGXcndK5Yz0ghNnMWBaEXwdWzsLKBf1bLqqfAnYC7HbZTxyGnhRQ2yb5c7bjrLFdo43LOOAo31iWjxP3aXcnVI5Y884IDiDqhSCgoAQFAQEAAAAAAAAyIx/BBONpqn84oIAAAAASUVORK5CYII=", "text/plain": [ "# \"Q-Q plot of the residuals\", :xlabel => \"Normal theoretical quantiles\", :ylabel => \"Observed quantiles\", :term => [\"png\"]], @datasets=Hamster::Vector[#, @options=Hamster::Hash[:pointtype => 6, :notitle => true, :with => \"points\"]>, # true, :with => \"lines\"]>], @cmd=\"plot \">" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "require 'distribution'\n", "\n", "observed = model_fit.residuals.sort\n", "n = observed.length\n", "theoretical = (1..n).to_a.map { |t| Distribution::Normal.p_value(t.to_f/n.to_f) * model_fit.sigma}\n", "qq_plot = Plot.new([[theoretical, observed], with: 'points', pointtype: 6, notitle: true],\n", " ['x', with: 'lines', notitle: true],\n", " xlabel: 'Normal theoretical quantiles', ylabel: 'Observed quantiles',\n", " title: 'Q-Q plot of the residuals')\n", "qq_plot.term('png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In general, we see that the model residuals satisfy the normality assumption to a reasonable extent." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Estimation results\n", "\n", "Let's look at the estimated model parameters, and see what they reveal about the data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Fixed effects\n", "\n", "Looking at the fixed effects coefficients it is striking that the estimate corresponding to the effect of the blog post text length is almost zero. Possibly, the length of a blog post has practically no effect on how many comments the blog post will receive." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Obtained fixed effects coefficient estimates:\n", "{:intercept=>0.18944500002662124, :log_host_comments_avg=>0.7861228519075351, :host_trackbacks_avg=>0.05579925647682859, :length=>2.580329980620255e-05, :has_parent_with_comments_lvl_yes=>-0.48841082096511457}\n" ] } ], "source": [ "puts \"Obtained fixed effects coefficient estimates:\"\n", "puts model_fit.fix_ef" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The directionality of the obtained estimates implies that blog posts hosted on websites with a high average of comments per blog post, also tend to have more comments. Moreover, blog posts which have parent blog posts, tend to have fewer comments. \n", "The effects of the average number of trackbacks per blog post on the hosting website and the blog post length seem rather small in magnitude." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Based on the [Wald Z test statistics](https://en.wikipedia.org/wiki/Wald_test#Test_on_a_single_parameter), we can carry out hypotheses tests for each fixed effects terms $\\beta_{i}$, testing the null $H_{0} : \\beta_{i} = 0$ against the alternative $H_{a} : \\beta_{i} \\neq 0$.\n", "\n", "*Note:* The Wald methods for $p$-values and confidence intervals are not absolutely trustworthy and should be treated with caution, as pointed out in [this blog post](http://agisga.github.io/MixedModels_p_values_and_CI/).\n", "\n", "The corresponding (approximate) p-values are obtained with:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{:intercept=>0.0, :log_host_comments_avg=>0.0, :host_trackbacks_avg=>5.723421736547607e-12, :length=>0.0, :has_parent_with_comments_lvl_yes=>3.859734754030342e-08}" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model_fit.fix_ef_p(method: :wald)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Interestingly, all $p$-values are tiny, implying that each predictor has a very strong linear relationship with the response variable.\n", "\n", "However, we should be careful with our conclusions. The very small $p$-values can also be explained by the large sample size (>20000 observations)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also look at Wald confidence intervals for the fixed effects coefficient estimates, which are in general more informative than $p$-values." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47139052686000 rows: 5 cols: 3
lower95upper95coef
intercept0.151616976164225870.22727302388901660.18944500002662124
log_host_comments_avg0.77569048814320870.79655521567186140.7861228519075351
host_trackbacks_avg0.0399179624770390350.071680550476618140.05579925647682859
length2.1719996776498967e-052.988660283590613e-052.580329980620255e-05
has_parent_with_comments_lvl_yes-0.6625496425736548-0.31427199935657435-0.48841082096511457
" ], "text/plain": [ "\n", "#\n", " lower95 upper95 coef \n", " intercept 0.15161697 0.22727302 0.18944500 \n", "log_host_c 0.77569048 0.79655521 0.78612285 \n", "host_track 0.03991796 0.07168055 0.05579925 \n", " length 2.17199967 2.98866028 2.58032998 \n", "has_parent -0.6625496 -0.3142719 -0.4884108 \n" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "conf_int = model_fit.fix_ef_conf_int(level: 0.95, method: :wald)\n", "ci = Daru::DataFrame.rows(conf_int.values, order: [:lower95, :upper95], index: model_fit.fix_ef_names)\n", "ci[:coef] = model_fit.fix_ef.values\n", "ci" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We observe that none of the 95% confidence intervals contains 0, which suggest high statistical significance of the linear predictors (equivalent to $p$-values). Also, all of the intervals seem rather narrow." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Random effects\n", "\n", "We can look at the obtained random effects estimates (the values $b_d$ from the above equation, where $d\\in\\{mo, tu, we, th, fr, sa, su\\}$)." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Obtained random effects coefficient estimates:\n", "{:intercept_fr=>0.0, :intercept_mo=>0.0, :intercept_sa=>0.0, :intercept_su=>0.0, :intercept_th=>0.0, :intercept_tu=>0.0, :intercept_we=>0.0}\n" ] } ], "source": [ "puts \"Obtained random effects coefficient estimates:\"\n", "puts model_fit.ran_ef" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also look at the estimated correlation structure of the random effects:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
Daru::DataFrame:47139054735200 rows: 1 cols: 1
day
day0.0
" ], "text/plain": [ "\n", "#\n", " day \n", " day 0.0 \n" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model_fit.ran_ef_summary" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Interestingly, the estimates of the random effects coefficients and standard deviation are all zero!\n", "\n", "That is, we have a singular fit. Thus, our results imply that the variability between different days of the week is not large enough to justify non-zero random effects in this model. Practically, we can coclude that the day of the week on which a blog post is published has no effect on the number of comments that the blog post will receive in the first 24 hours." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusions\n", "\n", "This turned out to be an example of a degenerate model fit (random effects variance estimated to be zero), and we saw that `mixed_models` can handle degenerate fits very well. \n", "\n", "Of course we learned alot about the blog post data by doing this analysis. Here is a list of some of the findings.\n", "\n", "* Blog posts hosted on websites with a high average of comments per blog post, also tend to have more comments (obviously). \n", "\n", "* Blog posts which have parent blog posts, seem to have fewer comments (wonder why that is...).\n", "\n", "* The effect of the average number of trackbacks per blog post on the hosting website seems to have a very small effect on the number of comments of a given blog post.\n", "\n", "* The blog post text length has an extremely small positive effect on the number of comments.\n", "\n", "* All considered fixed effects predictor variables seem to have a significant influence on the number of comments that a blog post receives in the first 24 hours after publication (according to Wald Z tests).\n", "\n", "* The day of the week on which a blog post is published has practically no influence on the number of comments that the blog post will receive in the first 24 hours after publication." ] } ], "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 }