<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/nzcensus-on-github/"></div></div><div><a data-url="https://www.r-bloggers.com/nzcensus-on-github/" data-counturl="https://www.r-bloggers.com/nzcensus-on-github/" data-count="horizontal" data-text="nzcensus on GitHub" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://ellisp.github.io/blog/2016/08/04/nzcensus-gam-elastic-lm"> Peter's stats stuff - R</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Introduction</h2>
<p>A few months back the first, pre-CRAN versions of my <code>nzelect</code> package included some data from the New Zealand Census 2013. As noted in my <a target="_blank" href="http://ellisp.github.io/blog/2016/07/14/nzelect-cran">last post</a>, I’ve now split this into a separate <code>nzcensus</code> package, for ease of development and maintenance and to allow <code>nzelect</code> to fit within CRAN size restrictions.</p>
<p>The <code>nzcensus</code> package has a set of 60+ demographic variables aggregated at the level of meshblock, area unit, territorial authority and regional council. These variables have mostly been calculated by me from the counts in the Statistics New Zealand “meshblock” data set. For example, I’ve converted the count of individuals of European descent into a proportion of individuals who provided ethnicity information. This makes the <code>nzcensus</code> R package more analysis-ready than the published meshblock data because the calculated variables are suitable for comparisons across area groupings, which wasn’t the case with the straight counts.</p>
<p>The code that does these calculations is <a target="_blank" href="https://github.com/ellisp/nzelect">available on GitHub</a>. The preparation of <code>nzcensus</code> is still done as part of the <code>nzelect</code> project, although end users can install either or both separately. I’ve fixed a few errors in my census calculations found in the past few weeks but there’s always potential for more. In the event of any uncertainty, refer back to the definitive version on the Statistics New Zealand website.</p>
<h2>Interactive map of household income</h2>
<p>One of the variables that I didn’t have to calculate (because it is directly available in the original data) is median household income. Here’s an interactive map showing that variable by area unit. Zoom in and out to the map in the usual way for map interaction on your device; click on circles to get the name and income value for the area unit in which it is centred. Red means higher income; blue means lower. The darkest red colour means in the top 10% of area units’ by household income.</p>
<p><iframe height="430px" width="100%" src="http://ellisp.github.io/img/0046-leaflet-map.html"></iframe></p>
<p>Here’s the simple R code to create that map, drawing on the <a target="_blank" href="https://rstudio.github.io/leaflet/">leaflet R package</a>, which provides extremely easy to use access to the Leaflet JavaScript library:</p>
<figure>
<pre><code data-lang="r"><span># install the nzcensus package (note it is part of the nzelect GitHub repository):
</span><span>devtools</span><span>::</span><span>install_github</span><span>(</span><span>&quot;ellisp/nzelect/pkg2&quot;</span><span>)</span><span>
</span><span># load up packages:
</span><span>library</span><span>(</span><span>leaflet</span><span>)</span><span>
</span><span>library</span><span>(</span><span>nzcensus</span><span>)</span><span>
</span><span># remove Chatham Islands </span><span>tmp</span><span> </span><span>&lt;-</span><span> </span><span>AreaUnits2013</span><span>[</span><span>AreaUnits2013</span><span>$</span><span>WGS84Longitude</span><span> </span><span>></span><span> </span><span>0</span><span> </span><span>&amp;</span><span> </span><span>!</span><span>is.na</span><span>(</span><span>AreaUnits2013</span><span>$</span><span>MedianIncome2013</span><span>),</span><span> </span><span>]</span><span>
</span><span># create colour palette function
</span><span>pal</span><span> </span><span>&lt;-</span><span> </span><span>colorQuantile</span><span>(</span><span>&quot;RdBu&quot;</span><span>,</span><span> </span><span>NULL</span><span>,</span><span> </span><span>n</span><span> </span><span>=</span><span> </span><span>10</span><span>)</span><span>
</span><span># create labels for popups
</span><span>labs</span><span> </span><span>&lt;-</span><span> </span><span>paste0</span><span>(</span><span>tmp</span><span>$</span><span>AU_NAM</span><span>,</span><span> </span><span>&quot; $&quot;</span><span>,</span><span> </span><span>format</span><span>(</span><span>tmp</span><span>$</span><span>MedianIncome2013</span><span>,</span><span> </span><span>big.mark</span><span> </span><span>=</span><span> </span><span>&quot;,&quot;</span><span>))</span><span>
</span><span># draw map:
</span><span>leaflet</span><span>()</span><span> </span><span>%>%</span><span>
</span><span>addProviderTiles</span><span>(</span><span>&quot;CartoDB.Positron&quot;</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>addCircles</span><span>(</span><span>lng</span><span> </span><span>=</span><span> </span><span>tmp</span><span>$</span><span>WGS84Longitude</span><span>,</span><span> </span><span>lat</span><span> </span><span>=</span><span> </span><span>tmp</span><span>$</span><span>WGS84Latitude</span><span>,</span><span>
</span><span>color</span><span> </span><span>=</span><span> </span><span>pal</span><span>(</span><span>-</span><span>tmp</span><span>$</span><span>MedianIncome2013</span><span>),</span><span>
</span><span>popup</span><span> </span><span>=</span><span> </span><span>labs</span><span>,</span><span>
</span><span>radius</span><span> </span><span>=</span><span> </span><span>500</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>addLegend</span><span>(</span><span>
</span><span>pal</span><span> </span><span>=</span><span> </span><span>pal</span><span>,</span><span>
</span><span>values</span><span> </span><span>=</span><span> </span><span>-</span><span>tmp</span><span>$</span><span>MedianIncome2013</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Quantile of median&lt;br>household income&quot;</span><span>,</span><span>
</span><span>position</span><span> </span><span>=</span><span> </span><span>&quot;topleft&quot;</span><span>,</span><span>
</span><span>bins</span><span> </span><span>=</span><span> </span><span>5</span><span>)</span></code></pre>
</figure>
<h2>Modelling household income on demographic variables</h2>
<p>In a <a target="_blank" href="http://ellisp.github.io/blog/2016/06/05/bootstrap-cv-strategies">previous post</a> which became one of my most read posts, I had used an earlier version of this dataset to try out different validation methods for regression modelling strategies. I mentioned at the time that the modelling strategies I was using (full model, stepwise selection, dropping variance-inflating collinear variables) were less important for that particular purpose than the three different validation methods, but I think most of the take-up in comments and on Twitter showed interest in the modelling strategies.</p>
<p>To take this idea to its next step, I start edging towards what I think is a better modelling strategy for that data. In today’s post I start with some techniques which I think both are needed for the optimal end approach, but I don’t yet bring them together. The methods I’m using here are:</p>
<ul>
<li>Adding latitude and longitude as explanatory variables to adjust for the fact that spatially located observations are not independent of eachother, but spatially close variables are likely to have correlated random elements;</li>
<li>Lasso, ridge regression, and elastic net regularisation (a generalisation of the lasso and ridge regression) to deal with the multicollinearity in a better way than simply dropping some of the variables;</li>
<li>A generalized additive model which lets me take into account non-linear shaped relationships between some of the explanatory variables and the response variable.</li>
</ul>
<p>My ultimate approach to this kind of regression modelling for this data – which has to wait for a later post – will combine the above, plus imputation of missing values for some explanatory variables, using the <code>gamsel</code> approach.</p>
<p>I’m going to fit three different models:</p>
<ul>
<li>a linear model with all explanatory variables, fit using ordinary least squares;</li>
<li>a linear model with all explanatory variables, fit via elastic net regularisation which penalises for extra variables included in the model and shrinks estimated coefficients towards zero to compensate;</li>
<li>a generalised additive model that lets the explanatory variables with the stronger relationship to income to have curved, non-linear relationships that better fit the actual data; and allows the nuisance spatial relationship with income to be a smooth, flexible spline.</li>
</ul>
<h3>Setup</h3>
<p>First I load up the packages with the R functionality I’m going to need and get the data into convenient shape for this particular modelling exercise.</p>
<figure>
<pre><code data-lang="r"><span>#===================setup=======================
</span><span>library</span><span>(</span><span>ggplot2</span><span>)</span><span>
</span><span>library</span><span>(</span><span>scales</span><span>)</span><span>
</span><span>library</span><span>(</span><span>boot</span><span>)</span><span>
</span><span>library</span><span>(</span><span>dplyr</span><span>)</span><span>
</span><span>library</span><span>(</span><span>Hmisc</span><span>)</span><span> </span><span># for spearman2
</span><span>library</span><span>(</span><span>mgcv</span><span>)</span><span> </span><span># for gam
</span><span>library</span><span>(</span><span>caret</span><span>)</span><span> </span><span># for RMSE
</span><span>library</span><span>(</span><span>grid</span><span>)</span><span>
</span><span>library</span><span>(</span><span>stringr</span><span>)</span><span>
</span><span>library</span><span>(</span><span>ggrepel</span><span>)</span><span>
</span><span>library</span><span>(</span><span>glmnet</span><span>)</span><span>
</span><span>library</span><span>(</span><span>maps</span><span>)</span><span> </span><span># for country borders on the smoothing plot
</span><span>
</span><span>set.seed</span><span>(</span><span>123</span><span>)</span><span>
</span><span># install nzcensus if not already done
# devtools::install_github(&quot;ellisp/nzelect/pkg2&quot;)
</span><span>library</span><span>(</span><span>nzcensus</span><span>)</span><span>
</span><span>au</span><span> </span><span>&lt;-</span><span> </span><span>AreaUnits2013</span><span> </span><span>%>%</span><span>
</span><span># drop the columns with areas' code and name, and the redundant coordinate system
</span><span> </span><span>select</span><span>(</span><span>-</span><span>AU2014</span><span>,</span><span> </span><span>-</span><span>AU_NAM</span><span>,</span><span> </span><span>-</span><span>NZTM2000Easting</span><span>,</span><span> </span><span>-</span><span>NZTM2000Northing</span><span>)</span><span> </span><span>%>%</span><span>
</span><span># drop some intrinsically collinear variables that would be exactly collinear
</span><span> </span><span># if it weren't for rounding error and confidentialisation:
</span><span> </span><span>select</span><span>(</span><span>-</span><span>PropWorked40_49hours2013</span><span>,</span><span> </span><span>-</span><span>Prop35to39_2013</span><span>,</span><span> </span><span>-</span><span>PropFemale2013</span><span>)</span><span>
</span><span># give meaningful rownames, helpful for some diagnostic plots later
</span><span>row.names</span><span>(</span><span>au</span><span>)</span><span> </span><span>&lt;-</span><span> </span><span>AreaUnits2013</span><span>$</span><span>AU_NAM</span><span>
</span><span># remove some repetition from the variable names
</span><span>names</span><span>(</span><span>au</span><span>)</span><span> </span><span>&lt;-</span><span> </span><span>gsub</span><span>(</span><span>&quot;_2013&quot;</span><span>,</span><span> </span><span>&quot;&quot;</span><span>,</span><span> </span><span>names</span><span>(</span><span>au</span><span>))</span><span>
</span><span>names</span><span>(</span><span>au</span><span>)</span><span> </span><span>&lt;-</span><span> </span><span>gsub</span><span>(</span><span>&quot;2013&quot;</span><span>,</span><span> </span><span>&quot;&quot;</span><span>,</span><span> </span><span>names</span><span>(</span><span>au</span><span>))</span><span>
</span><span>names</span><span>(</span><span>au</span><span>)</span><span> </span><span>&lt;-</span><span> </span><span>gsub</span><span>(</span><span>&quot;Prop&quot;</span><span>,</span><span> </span><span>&quot;&quot;</span><span>,</span><span> </span><span>names</span><span>(</span><span>au</span><span>))</span><span>
</span><span># restrict to areas with no missing data. An improvement for later
# is to use imputation and include this step within the validation.
</span><span>au</span><span> </span><span>&lt;-</span><span> </span><span>au</span><span>[</span><span>complete.cases</span><span>(</span><span>au</span><span>),</span><span> </span><span>]</span><span>
</span><span># give the data a generic name for ease of copying and pasting
</span><span>the_data</span><span> </span><span>&lt;-</span><span> </span><span>au</span><span>
</span><span># we need a dummy variable for the Chathams because it's extreme value of longitude
# makes any spatial variables otherwise highly problematic.
# the_data$chathams &lt;- the_data$WGS84Longitude &lt; 100
</span><span>the_data</span><span> </span><span>&lt;-</span><span> </span><span>the_data</span><span>[</span><span>the_data</span><span>$</span><span>WGS84Longitude</span><span> </span><span>></span><span> </span><span>100</span><span>,</span><span> </span><span>]</span><span>
</span><span>names</span><span>(</span><span>the_data</span><span>)</span><span> </span><span>&lt;-</span><span> </span><span>make.names</span><span>(</span><span>names</span><span>(</span><span>the_data</span><span>))</span></code></pre>
</figure>
<h3>Linear model</h3>
<p>My first and simplest method is one of those I used in the previous post: linear regression, estimated with ordinary least squares, using all the explanatory variables.</p>
<p>In comparing these approaches I am always going to use the same validation approach. All three validation methods gave similar results, so I’m going to pick one that’s easy to generalise, the so-called “simple bootstrap” validation method. In this approach, a re-sample of data is taken from the original data with replacement, the same number of rows as originally (hence some points sampled twice, some not at all). The model is fit based on this re-sample, and that model fit is used to predict the response variable for the entire original dataset (not just the “out of bag” samples – that’s a different approach). I collect the average root mean square error of those predictions from doing this resampling 99 times.</p>
<p>Because this classic model forms a reference point for alternative methods, it’s important to consider the classic diagnostic graphics too. Here they are:</p>
<p><img alt="diagnostics" src="http://ellisp.github.io/img/0046-lm-diagnostics.svg"></p>
<p>There’s no problems with outliers or leverage (noting that in the earlier preparation code I’d excluded the Chatham Islands from the dataset to make the spatial aspect of the problem tractable). There’s a bit of curvature in the residuals, indicating a possible structural problem that later on we will use the generalized additive model to address. There’s also some non-normality in the residuals – the QQ plot (top right of the four) shows that towards the upper and lower extremes of the residuals’ distribution the actual values of standardized residuals are larger in absolute magnitude than the theoretical values they should have from a Normal distribution. This isn’t disastrous, but it does mean we need to take care in inference, and use methods such as the bootstrap for all of our inference as methods that rely on distributional assuptions will be suspect.</p>
<p>Here’s the code for fitting this simple model. As will be the case for all three of our approaches, we start by defining a function that</p>
<ul>
<li>fits the model to a resample with replacement of the data, defined by the argument <code>i</code>;</li>
<li>uses that model to predict the values of income from the full original data;</li>
<li>returns the root mean square errors of those predictions.</li>
</ul>
<p>This function can then be used by the <code>boot</code> function to perform the bootstrap validation.</p>
<figure>
<pre><code data-lang="r"><span>fit_lm</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>data</span><span>,</span><span> </span><span>i</span><span>){</span><span>
</span><span>mod1</span><span> </span><span>&lt;-</span><span> </span><span>lm</span><span>(</span><span>MedianIncome</span><span> </span><span>~</span><span> </span><span>.</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>data</span><span>[</span><span>i</span><span>,</span><span> </span><span>])</span><span>
</span><span># use the model based on resample on the original data to estimate how
</span><span> </span><span># good or not it is:
</span><span> </span><span>RMSE</span><span>(</span><span>predict</span><span>(</span><span>mod1</span><span>,</span><span> </span><span>newdata</span><span> </span><span>=</span><span> </span><span>data</span><span>),</span><span> </span><span>data</span><span>$</span><span>MedianIncome</span><span>)</span><span>
</span><span>}</span><span>
</span><span># use bootstrap validation for an unbiaased estimate of root mean square error
</span><span>rmses_lm_boot</span><span> </span><span>&lt;-</span><span> </span><span>boot</span><span>(</span><span>the_data</span><span>,</span><span> </span><span>statistic</span><span> </span><span>=</span><span> </span><span>fit_lm</span><span>,</span><span> </span><span>R</span><span> </span><span>=</span><span> </span><span>99</span><span>)</span><span>
</span><span>lm_rmse</span><span> </span><span>&lt;-</span><span> </span><span>mean</span><span>(</span><span>rmses_lm_boot</span><span>$</span><span>t</span><span>)</span><span>
</span><span># save a single version of this model for later
</span><span>mod_lm</span><span> </span><span>&lt;-</span><span> </span><span>lm</span><span>(</span><span>MedianIncome</span><span> </span><span>~</span><span> </span><span>.</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>the_data</span><span>)</span><span>
</span><span># Check out classic diagnostic plots:
</span><span>par</span><span>(</span><span>mfrow</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>2</span><span>,</span><span> </span><span>2</span><span>))</span><span>
</span><span>plot</span><span>(</span><span>mod_lm</span><span>)</span></code></pre>
</figure>
<p>The average bootstrapped root mean square error from this method is 2015.</p>
<h3>Elastic net regularization</h3>
<p><a target="_blank" href="https://en.wikipedia.org/wiki/Elastic_net_regularization">Elastic net regularization</a> is a generalisation of the “ridge regression” and “lasso” methods of estimating parameters of a linear model. Both ridge regression and lasso estimation provide extra penalties for the magnitude of coefficients in the linear predictor; the “lasso” penalises based on the sum of the absolute value of estimated coefficients, and the “ridge” penalises based on the sum of squared values. In both cases the net effect is to shrink estimated coefficients towards zero. It’s a far superior method for model building than dropping variables based on their variable inflation factors, or stepwise regression, the two methods I used for demonstration purposes in my previous posts. In particular, all sorts of problems for subsequent inference, such as bias towards significance of the parameters remaining in the model, are reduced or eliminated altogether.</p>
<p>The elastic net approach is defined by two hyperparameters, <code>alpha</code> and <code>lambda</code>. When <code>alpha = 1</code> the elastic net is equivalent to the lasso; when <code>alpha = 0</code> it is equivalent to ridge regression. For any given value of <code>alpha</code>, the value of <code>lambda</code> indicates how much penalty is given to the size of coefficients in the linear predictor. The <code>cv.glmnet</code> function uses cross validation to give the best value of <code>lambda</code> for any given <code>alpha</code>; a bit of extra work is needed to simultaneously choose alpha. In the code below I use <code>cv.glmnet</code> with nine different values of <code>alpha</code>:</p>
<figure>
<pre><code data-lang="r"><span># A lasso, ridge rigression, or elastic net (which combines the two) is a way
# of dealing with the collinearity by forcing some coefficients to shrink (possibly to zero)
# while doing minimal damage to the inferential qualities of the rest and to the overall model fit.
</span><span>
</span><span># First we need to decide between ridge regression and lasso or elastic net (between the two)
# define folds for cross validation so can check the impact of different values of alpha </span><span>set.seed</span><span>(</span><span>124</span><span>)</span><span>
</span><span>foldid</span><span> </span><span>&lt;-</span><span> </span><span>sample</span><span>(</span><span>1</span><span>:</span><span>10</span><span>,</span><span> </span><span>nrow</span><span>(</span><span>the_data</span><span>),</span><span> </span><span>replace</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>)</span><span>
</span><span># separate out the explanatory from response variable
# for when using lasso and ridge regression
</span><span>X</span><span> </span><span>&lt;-</span><span> </span><span>the_data</span><span> </span><span>%>%</span><span> </span><span>select</span><span>(</span><span>-</span><span>MedianIncome</span><span>)</span><span>
</span><span>Y</span><span> </span><span>&lt;-</span><span> </span><span>the_data</span><span>$</span><span>MedianIncome</span><span> </span><span>cv_results</span><span> </span><span>&lt;-</span><span> </span><span>data_frame</span><span>(</span><span>lambda</span><span> </span><span>=</span><span> </span><span>numeric</span><span>(),</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>numeric</span><span>(),</span><span> </span><span>mcve</span><span> </span><span>=</span><span> </span><span>numeric</span><span>())</span><span>
</span><span>alphas</span><span> </span><span>&lt;-</span><span> </span><span>seq</span><span>(</span><span>from</span><span> </span><span>=</span><span> </span><span>0</span><span>,</span><span> </span><span>to</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>length.out</span><span> </span><span>=</span><span> </span><span>9</span><span>)</span><span>
</span><span>for</span><span>(</span><span>i</span><span> </span><span>in</span><span> </span><span>alphas</span><span>){</span><span>
</span><span>cvfit</span><span> </span><span>&lt;-</span><span> </span><span>cv.glmnet</span><span>(</span><span>as.matrix</span><span>(</span><span>X</span><span>),</span><span> </span><span>Y</span><span>,</span><span> </span><span>foldid</span><span> </span><span>=</span><span> </span><span>foldid</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>i</span><span>)</span><span>
</span><span>tmp</span><span> </span><span>&lt;-</span><span> </span><span>data_frame</span><span>(</span><span>lambda</span><span> </span><span>=</span><span> </span><span>cvfit</span><span>$</span><span>lambda</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>i</span><span>,</span><span> </span><span>mcve</span><span> </span><span>=</span><span> </span><span>cvfit</span><span>$</span><span>cvm</span><span>)</span><span> </span><span>cv_results</span><span> </span><span>&lt;-</span><span> </span><span>rbind</span><span>(</span><span>cv_results</span><span>,</span><span> </span><span>tmp</span><span>)</span><span>
</span><span>}</span><span>
</span><span>arrange</span><span>(</span><span>cv_results</span><span>,</span><span> </span><span>mcve</span><span>)</span><span> </span><span># best alpha with this see is 0.75 but right combination of alpha and lambda
# works pretty well for any alpha
</span><span>
</span><span># For the graphic I take square root of mcve so it is back on same scale as RMSE used elsewhere in this post
</span><span> </span><span>ggplot</span><span>(</span><span>cv_results</span><span>,</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>lambda</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>sqrt</span><span>(</span><span>mcve</span><span>),</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>as.factor</span><span>(</span><span>round</span><span>(</span><span>alpha</span><span>,</span><span> </span><span>3</span><span>))))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>2</span><span>)</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>0.2</span><span>,</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>&quot;grey50&quot;</span><span>,</span><span> </span><span>aes</span><span>(</span><span>group</span><span> </span><span>=</span><span> </span><span>as.factor</span><span>(</span><span>round</span><span>(</span><span>alpha</span><span>,</span><span> </span><span>3</span><span>))))</span><span> </span><span>+</span><span>
</span><span>scale_x_log10</span><span>()</span><span> </span><span>+</span><span>
</span><span>coord_cartesian</span><span>(</span><span>ylim</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>1800</span><span>,</span><span> </span><span>4000</span><span>),</span><span> </span><span>xlim</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>5</span><span>,</span><span> </span><span>1000</span><span>))</span><span> </span><span>+</span><span>
</span><span>scale_colour_brewer</span><span>(</span><span>&quot;alpha&quot;</span><span>,</span><span> </span><span>palette</span><span> </span><span>=</span><span> </span><span>&quot;Greens&quot;</span><span>,</span><span> </span><span>guide</span><span> </span><span>=</span><span> </span><span>guide_legend</span><span>(</span><span>reverse</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>))</span><span> </span><span>+</span><span>
</span><span>ggtitle</span><span>(</span><span>&quot;Cross-validation to select hyper parameters\nin elastic net regression&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>scale_y_continuous</span><span>(</span><span>&quot;Square root of mean cross validation error&quot;</span><span>,</span><span> </span><span>label</span><span> </span><span>=</span><span> </span><span>comma</span><span>)</span><span> </span><span>+</span><span>
</span><span>theme</span><span>(</span><span>legend.position</span><span> </span><span>=</span><span> </span><span>&quot;right&quot;</span><span>)</span></code></pre>
</figure>
<p>Here’s the results in tabular form. Many combinations of alpha and lambda work well. Basically, there’s not much to choose from between ridge regression and lasso with this particular set of data, which has a nice moderate size of observations for a modest number of variables.</p>
<figure>
<pre><code data-lang="r"><span>></span><span> </span><span>arrange</span><span>(</span><span>cv_results</span><span>,</span><span> </span><span>mcve</span><span>)</span><span>
</span><span>Source</span><span>:</span><span> </span><span>local</span><span> </span><span>data</span><span> </span><span>frame</span><span> </span><span>[</span><span>740</span><span> </span><span>x</span><span> </span><span>3</span><span>]</span><span>
</span><span>lambda</span><span> </span><span>alpha</span><span> </span><span>mcve</span><span>
</span><span>&lt;</span><span>dbl</span><span>></span><span> </span><span>&lt;</span><span>dbl</span><span>></span><span> </span><span>&lt;</span><span>dbl</span><span>></span><span>
</span><span>1</span><span> </span><span>6.029743</span><span> </span><span>0.750</span><span> </span><span>4290529</span><span>
</span><span>2</span><span> </span><span>5.168351</span><span> </span><span>0.875</span><span> </span><span>4291482</span><span>
</span><span>3</span><span> </span><span>4.963227</span><span> </span><span>1.000</span><span> </span><span>4291871</span><span>
</span><span>4</span><span> </span><span>7.235692</span><span> </span><span>0.625</span><span> </span><span>4291920</span><span>
</span><span>5</span><span> </span><span>7.941163</span><span> </span><span>0.625</span><span> </span><span>4292204</span><span>
</span><span>6</span><span> </span><span>6.592893</span><span> </span><span>0.625</span><span> </span><span>4292411</span><span>
</span><span>7</span><span> </span><span>4.290856</span><span> </span><span>0.875</span><span> </span><span>4292494</span><span>
</span><span>8</span><span> </span><span>4.522307</span><span> </span><span>1.000</span><span> </span><span>4292497</span><span>
</span><span>9</span><span> </span><span>6.617636</span><span> </span><span>0.750</span><span> </span><span>4292529</span><span>
</span><span>10</span><span> </span><span>5.672259</span><span> </span><span>0.875</span><span> </span><span>4292609</span><span>
</span><span>..</span><span> </span><span>...</span><span> </span><span>...</span><span> </span><span>...</span></code></pre>
</figure>
<p>Here’s the graphic presentation of those results, making it clearer that for an appropriate value of <code>lambda</code>, pretty much any alpha gives an ok result.</p>
<p><img alt="alpha-lambda" src="http://ellisp.github.io/img/0046-alpha-lambda.svg"></p>
<p>Choosing <code>0.750</code> as the value of <code>alpha</code> and using cross-validation for each resampled dataset to choose <code>lambda</code>, we can now use bootstrap validation to check out the results of this modelling approach.</p>
<figure>
<pre><code data-lang="r"><span>fit_elastic</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>data</span><span>,</span><span> </span><span>i</span><span>){</span><span>
</span><span># i = sample(1:nrow(data), nrow(data), replace = TRUE)
</span><span> </span><span>Y_orig</span><span> </span><span>&lt;-</span><span> </span><span>data</span><span>$</span><span>MedianIncome</span><span>
</span><span>X_orig</span><span> </span><span>&lt;-</span><span> </span><span>as.matrix</span><span>(</span><span>select</span><span>(</span><span>data</span><span>,</span><span> </span><span>-</span><span>MedianIncome</span><span>))</span><span>
</span><span>data2</span><span> </span><span>&lt;-</span><span> </span><span>data</span><span>[</span><span>i</span><span>,</span><span> </span><span>]</span><span>
</span><span>Y_new</span><span> </span><span>&lt;-</span><span> </span><span>data2</span><span>$</span><span>MedianIncome</span><span>
</span><span>X_new</span><span> </span><span>&lt;-</span><span> </span><span>as.matrix</span><span>(</span><span>select</span><span>(</span><span>data2</span><span>,</span><span> </span><span>-</span><span>MedianIncome</span><span>))</span><span>
</span><span>lambda</span><span> </span><span>&lt;-</span><span> </span><span>cv.glmnet</span><span> </span><span>(</span><span>X_new</span><span>,</span><span> </span><span>Y_new</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.85</span><span>)</span><span>$</span><span>lambda.min</span><span>
</span><span>mod1</span><span> </span><span>&lt;-</span><span> </span><span>glmnet</span><span>(</span><span>X_new</span><span>,</span><span> </span><span>Y_new</span><span>,</span><span> </span><span>lambda</span><span> </span><span>=</span><span> </span><span>lambda</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.85</span><span>)</span><span>
</span><span># use the model based on resample on the original data to estimate how
</span><span> </span><span># good or not it is:
</span><span> </span><span>rmse</span><span> </span><span>&lt;-</span><span> </span><span>RMSE</span><span>(</span><span>predict</span><span>(</span><span>mod1</span><span>,</span><span> </span><span>newx</span><span> </span><span>=</span><span> </span><span>X_orig</span><span>),</span><span> </span><span>Y</span><span>)</span><span> </span><span>return</span><span>(</span><span>rmse</span><span>)</span><span>
</span><span>}</span><span>
</span><span>rmses_elastic_boot</span><span> </span><span>&lt;-</span><span> </span><span>boot</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>the_data</span><span>,</span><span> </span><span>statistic</span><span> </span><span>=</span><span> </span><span>fit_elastic</span><span>,</span><span> </span><span>R</span><span> </span><span>=</span><span> </span><span>99</span><span>)</span><span> </span><span># takes a few minutes
</span><span>elastic_rmse</span><span> </span><span>&lt;-</span><span> </span><span>mean</span><span>(</span><span>rmses_elastic_boot</span><span>$</span><span>t</span><span>)</span></code></pre>
</figure>
<p>The average bootstrapped root mean square error from this method is again 2015, showing that the improved estimates of individual coefficients has neither worsened nor improved the overall fit.</p>
<p>The characteristic of shrinkage methods is that the estimates of coefficients in the linear predictor are less in absolute size than are the ordinary least squares estimators. The relatively mild nature of this phenomenon in this particular case is shown in the graphic below. The straight diagonal line shows where the shrunk estimate of the coefficient from the elastic net equals that from ordinary least squares; when points are to the right of this line if positive, or to the left of it if negative, the elastic net regularized estimate is closer to zero than its ordinary least squares equivalent. The shrinkage is pretty mild.</p>
<p><img alt="coefs" src="http://ellisp.github.io/img/0046-compare-coefs.svg"></p>
<h3>Generalized additive model</h3>
<p>The third and final strategy to be tried today is a different structural model, allowing non-linear relations between the explanatory variables and the response (median household income) they are trying to predict. It turns out this has a fairly material impact on overall model fit. I’ll start with a plot showing those non-linear relations for the fifteen demographic variables most closely related to household income and the two spatial (latitude and longitude) variables:</p>
<p><img alt="curves" src="http://ellisp.github.io/img/0046-gam-relations.svg"></p>
<p>We can see a bit of curve in some of these relationships. For example, the relationship between the proportion of people who are partnered and income is decidedly non-linear, with an upwards surge in income for a small number of low-partnership areas that belies the general positive correlation. The relationship between the proportion of people in full time employment and income shows a definite S shape. The relationship between the proportion of people with Bachelors degrees and income gets steeper as the proportion gets higher; and so on.</p>
<p>The latitude and longitude plot at the bottom right of the image above is the most complex non-linear relationship and the plot is a little difficult to see, so here’s an enhanced version of it:</p>
<p><img alt="map" src="https://i1.wp.com/ellisp.github.io/img/0046-gam-spatial-residuals.png?w=456" data-recalc-dims="1"></p>
<p>This map is showing the relationship between location and income, after all the demographic variables are accounted for. Including a flexible spline for an interaction between longitude and latitude in a Generalized Additive Model like this is a good way of controlling for the nuisance factor of spatial autocorrelation, which if not accounted for is likely to lead to false inference, with the strength of evidence for relations between the data being overestimated (apparent relations can be just artefacts of how the spatial grouping was made).</p>
<p>In specifying models of this sort I often follow Frank Harrell’s informal advice of first determining how many degrees of freedom are available (between 1/20 and 1/10 the number of observations). Those degrees of freedom are then allocated to the various explanatory variables in the form of splines, with more degrees of freedom meaning a flexible curved relationship with income is possible and a single degree of freedom meaning a straight linear relationship. The explanatory variables with the strongest non-parametric relationship with the response variable get more degrees of freedoms. This approach works only after the analyst has made a commitment to include all examined variables in the model; letting the strength of relationship with the response determine whether an explanatory variable gets in at all is a form of stepwise regression and subject to all the problems that come from that.</p>
<p>To estimate the strength of the non parameteric relationships I use Harrell’s <code>spearman2</code> which regresses Y on the rank of X and the square of the rank of X. This powerful method allows non-monotonic relationships between variables. In the case of today’s data, it shows the strongest relationship to household income as being the proportion of people in full time employment, households with access to the internet, and so on as seen in the output below:</p>
<figure>
<pre><code data-lang="r"><span>></span><span> </span><span>sp</span><span> </span><span>&lt;-</span><span> </span><span>spearman2</span><span>(</span><span>MedianIncome</span><span> </span><span>~</span><span> </span><span>.</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>the_data</span><span>)</span><span>
</span><span>></span><span> </span><span>sp</span><span>[</span><span>order</span><span>(</span><span>-</span><span>sp</span><span>[</span><span> </span><span>,</span><span>6</span><span>])[</span><span>1</span><span>:</span><span>15</span><span>],</span><span> </span><span>]</span><span>
</span><span>rho2</span><span> </span><span>F</span><span> </span><span>df1</span><span> </span><span>df2</span><span> </span><span>P</span><span> </span><span>Adjusted</span><span> </span><span>rho2</span><span> </span><span>n</span><span>
</span><span>FullTimeEmployed</span><span> </span><span>0.7134528</span><span> </span><span>4436.8712</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.7132920</span><span> </span><span>1784</span><span>
</span><span>InternetHH</span><span> </span><span>0.5844727</span><span> </span><span>2506.5272</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.5842396</span><span> </span><span>1784</span><span>
</span><span>NoQualification</span><span> </span><span>0.4477222</span><span> </span><span>1444.6370</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.4474123</span><span> </span><span>1784</span><span>
</span><span>UnemploymentBenefit</span><span> </span><span>0.4309570</span><span> </span><span>1349.5736</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.4306377</span><span> </span><span>1784</span><span>
</span><span>Smoker</span><span> </span><span>0.4103768</span><span> </span><span>1240.2690</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.4100459</span><span> </span><span>1784</span><span>
</span><span>Partnered</span><span> </span><span>0.3866804</span><span> </span><span>1123.5000</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.3863363</span><span> </span><span>1784</span><span>
</span><span>Managers</span><span> </span><span>0.3854615</span><span> </span><span>1117.7367</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.3851166</span><span> </span><span>1784</span><span>
</span><span>Bachelor</span><span> </span><span>0.3729497</span><span> </span><span>1059.8773</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.3725978</span><span> </span><span>1784</span><span>
</span><span>SelfEmployed</span><span> </span><span>0.3666249</span><span> </span><span>1031.4988</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.3662695</span><span> </span><span>1784</span><span>
</span><span>NoMotorVehicle</span><span> </span><span>0.3587586</span><span> </span><span>996.9846</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.3583987</span><span> </span><span>1784</span><span>
</span><span>Unemployed</span><span> </span><span>0.3572286</span><span> </span><span>990.3698</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.3568679</span><span> </span><span>1784</span><span>
</span><span>Labourers</span><span> </span><span>0.3387476</span><span> </span><span>912.8864</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.3383766</span><span> </span><span>1784</span><span>
</span><span>Worked50_59hours</span><span> </span><span>0.3311916</span><span> </span><span>882.4404</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.3308163</span><span> </span><span>1784</span><span>
</span><span>Separated</span><span> </span><span>0.3122393</span><span> </span><span>809.0175</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.3118533</span><span> </span><span>1784</span><span>
</span><span>Maori</span><span> </span><span>0.3022017</span><span> </span><span>771.7465</span><span> </span><span>1</span><span> </span><span>1782</span><span> </span><span>0</span><span> </span><span>0.3018101</span><span> </span><span>1784</span></code></pre>
</figure>
<p>I then allocate flexible splines to the first 15 or so of those variables. The 4 x 4 grid of non-linear plots and the map of spatial effects that began this section was produced with the code below:</p>
<figure>
<pre><code data-lang="r"><span># see http://stackoverflow.com/questions/30627642/issue-with-gam-function-in-r
# for why the data needs to be specified in both terms() and in gam()
</span><span>the_formula</span><span> </span><span>&lt;-</span><span> </span><span>terms</span><span>(</span><span>MedianIncome</span><span> </span><span>~</span><span> </span><span>s</span><span>(</span><span>FullTimeEmployed</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>6</span><span>)</span><span> </span><span>+</span><span> </span><span>s</span><span>(</span><span>InternetHH</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>6</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>NoQualification</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>5</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>UnemploymentBenefit</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>5</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>Smoker</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>5</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>Partnered</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>5</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>Managers</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>4</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>Bachelor</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>4</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>SelfEmployed</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>4</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>NoMotorVehicle</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>4</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>Unemployed</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>3</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>Labourers</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>3</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>Worked50_59hours</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>3</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>Separated</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>3</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>Maori</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>3</span><span>)</span><span> </span><span>+</span><span>
</span><span>s</span><span>(</span><span>WGS84Longitude</span><span>,</span><span> </span><span>WGS84Latitude</span><span>)</span><span> </span><span>+</span><span>
</span><span>.</span><span>,</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>the_data</span><span>)</span><span>
</span><span>gam_model</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>the_formula</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>the_data</span><span>)</span><span>
</span><span># grid of nonlinear relations
</span><span>par</span><span>(</span><span>bty</span><span> </span><span>=</span><span> </span><span>&quot;l&quot;</span><span>,</span><span> </span><span>mar</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>5</span><span>,</span><span>4</span><span>,</span><span> </span><span>2</span><span>,</span><span> </span><span>1</span><span>))</span><span>
</span><span>plot</span><span>(</span><span>gam_model</span><span>,</span><span> </span><span>residuals</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>pages</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>shade</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>seWithMean</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>ylab</span><span> </span><span>=</span><span> </span><span>&quot;&quot;</span><span>)</span><span>
</span><span>grid.text</span><span>(</span><span>&quot;Impact of area average variables on median income by area unit (New Zealand census 2013)&quot;</span><span>,</span><span> </span><span>0.5</span><span>,</span><span> </span><span>0.99</span><span>,</span><span>
</span><span>gp</span><span> </span><span>=</span><span> </span><span>gpar</span><span>(</span><span>fontfamily</span><span> </span><span>=</span><span> </span><span>&quot;myfont&quot;</span><span>))</span><span>
</span><span># map of spatial impacts	</span><span>par</span><span>(</span><span>bty</span><span> </span><span>=</span><span> </span><span>&quot;l&quot;</span><span>,</span><span> </span><span>family</span><span> </span><span>=</span><span> </span><span>&quot;myfont&quot;</span><span>,</span><span> </span><span>fg</span><span> </span><span>=</span><span> </span><span>&quot;white&quot;</span><span>)</span><span>
</span><span>plot</span><span>(</span><span>gam_model</span><span>,</span><span> </span><span>shade</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>select</span><span> </span><span>=</span><span> </span><span>16</span><span>,</span><span> </span><span>rug</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>se</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>,</span><span> </span><span>scheme</span><span> </span><span>=</span><span> </span><span>2</span><span>,</span><span> </span><span>col</span><span> </span><span>=</span><span> </span><span>topo.colors</span><span>(</span><span>100</span><span>),</span><span> </span><span>pch</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>ylab</span><span> </span><span>=</span><span> </span><span>&quot;Latitude&quot;</span><span>,</span><span> </span><span>xlab</span><span> </span><span>=</span><span> </span><span>&quot;Longitude&quot;</span><span>,</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;Spatial pattern in regression of income\non demographic area variables&quot;</span><span>)</span><span>
</span><span>map</span><span>(</span><span>add</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>col</span><span> </span><span>=</span><span> </span><span>&quot;grey75&quot;</span><span>)</span></code></pre>
</figure>
<p>The bootstrap validation of this method is straightforward and follows the same pattern as the two previous methods:</p>
<figure>
<pre><code data-lang="r"><span>fit_gam</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>data</span><span>,</span><span> </span><span>i</span><span>){</span><span>
</span><span>mod1</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>the_formula</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>data</span><span>[</span><span>i</span><span>,</span><span> </span><span>])</span><span>
</span><span># use the model based on resample on the original data to estimate how
</span><span> </span><span># good or not it is:
</span><span> </span><span>RMSE</span><span>(</span><span>predict</span><span>(</span><span>mod1</span><span>,</span><span> </span><span>newdata</span><span> </span><span>=</span><span> </span><span>data</span><span>),</span><span> </span><span>data</span><span>$</span><span>MedianIncome</span><span>)</span><span>
</span><span>}</span><span>
</span><span>rmses_gam_boot</span><span> </span><span>&lt;-</span><span> </span><span>boot</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>the_data</span><span>,</span><span> </span><span>statistic</span><span> </span><span>=</span><span> </span><span>fit_gam</span><span>,</span><span> </span><span>R</span><span> </span><span>=</span><span> </span><span>99</span><span>)</span><span> </span><span># takes a few minutes
</span><span>gam_rmse</span><span> </span><span>&lt;-</span><span> </span><span>mean</span><span>(</span><span>rmses_gam_boot</span><span>$</span><span>t</span><span>)</span></code></pre>
</figure>
<p>The root mean square error from this method is 1728, a very material improvement in model performance; much less than 2015.</p>
<p>The spatial spline (which was allowed to be as flexible as it needed) uses 27 effective degrees of freedom and the other splines another 30 (plus 60 or so from the parametric effects), taking this model to the most flexibility reasonable from the 1775 data points.</p>
<h3>Conclusions</h3>
<p>Here’s those overall results:</p>
<figure>
<pre><code data-lang="r"><span>></span><span> </span><span>data.frame</span><span>(</span><span>lm_rmse</span><span>,</span><span> </span><span>elastic_rmse</span><span>,</span><span> </span><span>gam_rmse</span><span>)</span><span>
</span><span>lm_rmse</span><span> </span><span>elastic_rmse</span><span> </span><span>gam_rmse</span><span>
</span><span>1</span><span> </span><span>2015.34</span><span> </span><span>2014.871</span><span> </span><span>1728.254</span></code></pre>
</figure>
<p>The elastic net approach gives us more confidence in the estimates of individual coefficients but in this case does not materially improve model performance. Allowing curvature in the relations improved the model but does not address the collinearity problems.</p>
<p>As foreshadowed earlier, the next improvements to make will be to combine these two methods – elastic net regularization for feature selection and shrinkage, and generalized additive models to allow non-linear relationships. The <code>gamsel</code> R package enables this and I’ll return to it in a later post; perhaps at the same time as addressing an improved method of dealing with missing data (imputation rather than dropping altogether).</p>
<p>For now, this serves as a sequel to my earlier piece on model validation, turning the attention now onto model strategies themselves; and to demonstrate a taster of the uses the <code>nzcensus</code> package can be turned to.</p>
<p>The <code>nzcensus</code> package is still experimental, use with caution. Please file bugs or enhancement requests as <a target="_blank" href="https://github.com/ellisp/nzelect/issues">issues on GitHub</a>.</p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="https://www.r-bloggers.com/nzcensus-on-github/"></div></div><div><a data-url="https://www.r-bloggers.com/nzcensus-on-github/" data-counturl="https://www.r-bloggers.com/nzcensus-on-github/" data-count="vertical" data-text="nzcensus on GitHub" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://ellisp.github.io/blog/2016/08/04/nzcensus-gam-elastic-lm"> Peter's stats stuff - R</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/7N3ejJrH5dE"></a></div>