pfh blog2017-11-08T05:34:59ZIdeas for free software projects, politics, philosophy, amateur psychology.Paul Harrisonpfh@logarithmic.netTopconfects talk2017-11-08T05:34:59Z2017-11-08T05:34:59Zhttp://www.logarithmic.net/pfh/blog/01510119299
<p>I gave an informal talk today about my <a href="http://logarithmic.net/topconfects/">Topconfects</a> R package. If you do RNA-seq Differential Expression analysis it may be of interest.
<p><a href="http://www.logarithmic.net/pfh-files/blog/01510119299/informal-talk.html">>> View slides</a>Scatter plots with density quartiles2017-10-28T03:55:40Z2017-10-28T03:55:40Zhttp://www.logarithmic.net/pfh/blog/01509162940
<p>I think this is a better way to show density in scatter plots.
<p><a href="http://www.logarithmic.net/pfh-files/blog/01509162940/scatter.html">>> Read about it.
<p><img src="http://www.logarithmic.net/pfh-files/blog/01509162940/scatter.png"></a>
<p>A Bayesian walks into a bar and observes that a statistical hypothesis has been rejected2017-06-24T06:44:35Z2017-06-24T06:44:35Zhttp://www.logarithmic.net/pfh/blog/01498286675
<p>A classical statistical test controls the probability P(R|H0) that H0 will be falsely rejected. Conventionally this is controlled to be α=0.05 or α=0.01.
<p>( As a small complication, H0 may be a set of hypotheses. α is the probability of rejection of the hypothesis in the set that is most likely to be rejected. A Bayesian may believe some of the hypotheses are more likely than others, but whatever they believe they will still have that P(R|H0) ≤ α. )
<p>A Bayesian, in order to update their belief in the odds of H1 as opposed to H0, P(H1)/P(H0), also needs to know the statistical power, P(R|H1), the probability that H0 will be rejected if it is actually false. The
<a href="https://en.wikipedia.org/wiki/Bayes_factor">Bayes Factor</a>
is then at least P(R|H1)/P(R|H0), and the odds of the competing hypotheses can be conservatively updated as:
<p>P(H1|R)/P(H0|R) = P(R|H1)/P(R|H0) * P(H1)/P(H0)
<p>So we have the curious conclusion that a Bayesian will pay more heed to a statistical test they believe to have high statistical power. If the Bayesian believes that H0 had little chance of being correctly rejected if false, they will be surprised by the rejection but it will not update their beliefs much.
<p>( The classical test sought only to reject H0, not to confirm H1. If the test is rejects H0, perhaps the Bayesian should consider that their H1 was not sufficiently broad an alternative hypothesis. )
<p><br>Confidence intervals are increasingly accepted as a superior alternative to p-values. A Bayesian argument for this can be given: A confidence interval gives an indication of the accuracy to which an effect has been measured, and hence the statistical power. A Bayesian may use a confidence interval to update their beliefs immediately, whereas they would require further information if only provided with a p-value. ( Leaving aside the technical distinction between confidence intervals and credible intervals, which does not exist for simple tests such as the t-test. )
<p><br>( The probabilities above are Bayesian -- personal and subjective quantities used to make decisions. If we were talking frequency probabilities, the above would earn a <a href="https://link.springer.com/article/10.1007%2Fs10654-016-0149-3">"No!"</a> )
<p><br>Melbourne Datathon 2017 - my Kaggle entry2017-05-29T23:28:08Z2017-05-29T23:28:08Zhttp://www.logarithmic.net/pfh/blog/01496100488
<p><ul><li><a href="http://www.logarithmic.net/pfh-files/blog/01496100488/make_db.R">make_db.R</a> - load data into SQLite</li></ul>
<ul><li><a href="http://www.logarithmic.net/pfh-files/blog/01496100488/predict.R">predict.R</a> - make predictions (except for Stan model)</li></ul>
<ul><li><a href="http://www.logarithmic.net/pfh-files/blog/01496100488/sparse_glm_csr.stan">sparse_glm_csr.stan</a> - Stan model code</li></ul>
<ul><li><a href="http://www.logarithmic.net/pfh-files/blog/01496100488/stan.R">stan.R</a> - make input for Stan model, read results and make final prediction</li></ul>
<p>This is a description of my
<a href="https://inclass.kaggle.com/c/dsm2017">Melbourne Datathon 2017 Kaggle</a>
entry, which came third. See also: <a href="https://inclass.kaggle.com/c/dsm2017/forums/t/33996/my-solution-to-datathon-2017">Yuan Li's winning entry</a>.
<p>The data this year was drug purchases at pharmacies from 2011-2015. The objective was then to maximize the
<a href="https://en.wikipedia.org/wiki/Receiver_operating_characteristic#Area_under_the_curve">ROC AUC</a>
when predicting the probability of a person purchasing diabetes medication in 2016.
<p>This entry spent quite a lot of time at the top of the leaderboard, however this was merely because the experienced Kaggle players only submitted in the last couple of days. On the plus side, I guess this gave new players a chance to compete with each other rather than be faced with a seemlingly insurmountable AUC from the start.
<p>One thing I wondered was how much I should discuss my entry with other people or keep it secret. Well it was much more fun when I chatted than when I keept silent, and I needn't have worried anyway. The whole thing was a great learning exercise, I've learned some things sure to be relevant to my own work.
<p><br><b>Features</b>
<p>All models I used followed a pattern of predicting the binary outcome for each patient (diabetes medication purchased in 2016) from a collection of predictor variables. These predictors were stored in a large, sparse matrix.
<p>I actually used two predictor matrices, one for logistic regression prediction, and a larger one for tree prediction with XGBoost. XGBoost seems able to handle a very large number of predictor variables.
<p>The logistic regression predictors were:
<p><ul><li>Purchased a diabetes drug prior to 2016.</li></ul>
<ul><li>Year of birth, missing values set to 1921.</li></ul>
<ul><li>Year of birth, missing values set to 2016.</li></ul>
<ul><li>Known to be female.</li></ul>
<ul><li>Known to be male.</li></ul>
<ul><li>Patient location code.</li></ul>
<ul><li>For each drug, number of purchases, transformed log(1+x).</li></ul>
<ul><li>For each drug, number of purchases in 2015, transformed log(1+x).</li></ul>
<ul><li>For each chronic illness, number of drug purchases, transformed log(1+x).</li></ul>
<ul><li>For each chronic illness, number of drug purchases in 2015, transformed log(1+x).</li></ul>
<ul><li>For each prescriber, 1 if had a prescription from prescriber else 0.</li></ul>
<ul><li>For each store, proportion of purchases from store (ie rows normalized to sum to 1).</li></ul>
<p>log(1+x) transformation worked better than either using the count as-is or an indivator variable of whether the count was non-zero. I view it as an intermediate between these extremes.
<p>Additionally the XGBoost predictor matrix contained number of drug purchases in the second half of 2015, from 2014, from 2013, and from 2012.
<p>All features were standardized to have standard deviation one.
<p><br><b>Principal Component augmentation</b>
<p>I used the R package
<a href="https://cran.r-project.org/web/packages/irlba/index.html">irlba</a>
to augment the feature matrices with their principal components. Irlba is able to perform SVD on large sparse matrices.
<p>I added 10 principal component columns to the matrix for tree models, and 20 to the matrix for logistic regression models. The main limitation here was that irlba started breaking or running a very long time when asked for more than this.
<p>One interesting thing here is that this adds an aspect of unsupervised learning that includes the test-set patients.
<p><br><b>Bootstrap aggregation (<a href="https://en.wikipedia.org/wiki/Bootstrap_aggregating">"bagging"</a>)</b>
<p>I wanted to take model uncertainty into account. If one possible model makes a very strong prediction, and another makes a weaker prediction, I wanted an intermediate prediction. The most likely model might not well characterize the full spread of possibilities.
<p>To this end, each model was bootstrapped 60 times. Probabilities from these 60 bootstraps were averaged. The idea here is to obtain marginal probabilities, integrating over the posterior distribution of models. (Someone more familiar with the theory of bootstrapping could maybe tell me if this is a correct use of it. It certainly improved my AUC.)
<p>Boostrapping is based on sampling with replacement. My bootstrapper function had one slightly fancy feature: overall, each data point is used the same number of times. Within this constraint, it mimics the distribution produced by sampling with replacement as closely as possible. I call this a "balanced bootstrap".
<p>Since bootstraps always leave out about a third of the data, averaging over predictions from the left-out data provides a way to get probabilities for training-set patients that are not over-fitted. These are then used when blending different models. My "balanced bootstrap" scheme ensured every training-set patient was left out of the same proportion of bootrstaps, and I could dial the whole thing down to 2 bootstraps and still have it work if I wanted.
<p><br><b>Blending models</b>
<p>The probabilities produced by different models were linearly blended so as to maximize AUC. I again used my bootstrapper to perform this, so I had a good prediction of how accurate it would be without needing to make a Kaggle submission. (Once you have a hammer...)
<p><br>On to the models themselves.
<p><br><b>Regularized logistic regression with glmnet</b>
<p><a href="https://cran.r-project.org/web/packages/glmnet/index.html">glmnet</a>
fits generalized linear models with a mixture of L2 (ridge) and L1 (lasso) regularization.
I used glmnet to perform logistic regression.
<p>glmnet's alpha parameter determines the mixture of L2 and L1 regression. alpha=0 means entirely L2 and alpha=1 means entirely L1. I found a tiny bit of L1 improved predictions slightly over purely L2 regularization. In the end I made one set of predictions using alpha=0.01 and one with alpha=0.
<p>A neat feature of glmnet is that it fits a range of regularization amounts (lambda parameter) in one go, producing a path through coefficient space. As described above, I am always fitting models to bootstrapped data sets. I choose an optimal value for lambda each time to maximize the AUC for the patients left out by the bootstrap.
<p><br>I want to mention one of the strengths of glmnet but not so relevant to a Kaggle competition: if you turn up alpha close or equal to 1, and turn up lambda high enough, it will produce a <i>sparse</i> set of coefficients with few enough non-zero coefficients that the model is <i>interpretable</i>. Here I've used a low value of alpha to get as accurate predictions as possible, but if I wanted an interpretable model I would do the opposite.
<p><br><br><b>Gradient boosted decision trees with XGBoost</b>
<p>Speaking of non-interpretable models, my next and individually most accurate set of predictions came from
<a href="https://cran.r-project.org/web/packages/xgboost/index.html">XGBoost</a>.
<p>XGBoost creates a stack of decision trees, each tree additively refining the predicted log odds ratio for each patient. By piling up trees like this, and rather like glmnet, we have a path from high regularization to low regularization. Again I used the patients left out of the bootstrap sample to decide at what point along this path to stop.
<p>I used a learning rate of eta=0.1, and allowed a maximum tree depth of 12. XGBoost also has lambda and alpha parameters controlling L1 and L2 regularization. However the meaning of these parameters is different to glmnet: lambda is L2 regularization and alpha is L1 regularization. I'm not sure of the units here. I made two sets of predictions, one with lambda=50, alpha=10, and one with lambda=0, alpha=40. L1 regularization can stop to tree splitting before the maximum depth of 12 is reached (I think).
<p>Limiting the learning rate and number of trees also constitutes a form of regularization, even without these explicit regularization parameters. I'm not 100% clear on how all these parameters interact.
<p><br><b>Split glmnet</b>
<p>The most important predictor of buying diabetes drugs in 2016 was buying diabetes drugs prior to 2016. I split the data into two parts on this basis, and fitted glmnet models to each.
<p><br><b>XGBoost on top of glmnet</b>
<p>I tried running XGBoost starting from half the log odds ratio predicted by glmnet, as an attempt to blend these two approaches more intimately.
<p><br><b>Logistic regression with t distribution priors on coefficients with Stan</b>
<p><a href="http://mc-stan.org/">Stan</a>
is a language for writing Bayesian models. A Stan model can then be compiled to a C++ program that performs Hamiltonian Monte Carlo sampling.
<p>This was a late addition. I didn't bootstrap this because it does its own sampling from the posterior, just blended it in 10%-90% with the blend of the predictions described above.
<p>L2 regulatization can be viewed as a prior belief that coefficients have a normal distribution, and L1 regularization can be viewed as a prior belief that coefficients have a Laplace distribution. However my favourite distribution is Student's t distribution, so I wanted to try this as well. The t distribution looks like a normal distribution, but with fat tails. The Laplace distribution has fatter tails than the normal distribution, but the t distribution's fat tails are fatter still.
<p>I constructed a very simple Stan model for logistic regression with t distributed coefficients. To my surprise this worked, and once I figured out how to use the sparse matrix multiplication function it ran within a bearable time. I ran-in the sampler over 100 iterations, then sampled 100 models from the posterior distribution and averaged the predicted probabilities from these models.
<p>I was impressed by Stan, and if I were doing this exercise over would have spent more time on it.
<p>So there you go, that's my Melbourne Datathon 2017 Kaggle entry in all it's <strike>horror</strike> glory.
<p><br><b>Post-competition notes:</b>
<p>I came third on the final leaderboard.
<p><a href="https://inclass.kaggle.com/c/dsm2017/forums/t/33996/my-solution-to-datathon-2017">Yuan Li has described their winning solution</a>, which is xgboost based with much better feature engineering than I've used.
<p>The competition allowed me to submit two final entries, and I submitted one with the Stan model contribution and one without. The one without turned out to be the best of the two, by a hair.
<p><a href="http://www.logarithmic.net/pfh-files/blog/01496100488/prediction.nb.html">A small attempt at interpretable modelling of the data, using L1 regularized logistic regression.</a>
<p>Better at predicting cesation of diabetes medication than starting it:
<p><img src="http://www.logarithmic.net/pfh-files/blog/01496100488/is_it_useful_small.png">Diagrams of classical statistical procedures2017-04-01T22:58:04Z2017-04-01T22:58:04Zhttp://www.logarithmic.net/pfh/blog/01491087484
<img src="http://www.logarithmic.net/pfh-files/blog/01491087484/acceptance.svg">
<p><br>No matter what the truth may be, classical statistical procedures only ever reject it with some small specified probability α. This all-possible-worlds counterfactualism requires a kind of thinking filled with 90-degree turns, and it does my head in. Never fear, however, for I have devised the diagrams above to make what is going on clear. While looking at these diagrams I experience the sensation of understanding. I hope to now share this sensation with you.
<p>We shall be concerned with estimating a single parameter, or "effect". There is a true effect, and we have an estimate of it, but we know the estimate is inaccurate and have a probabilistic model describing exactly how inaccurate it will be, P(estimate|true). Here we will assume it has added t-distributed noise with zero mean and known scale.
<p>Consider a diagram with the true value of a parameter as the x-axis, and the estimate as the y-axis. We will color in all points that are non-rejected.
<p>Such a diagram represents a <b>valid</b> procedure if, for each true value, the non-rejected region will contain the estimate with probability 1-α. So in assessing validity we look at each <b>vertical</b> slice of the diagram.
<p>To <b>apply</b> the diagram, given an estimated effect, we look at the corresponding <b>horizontal</b> slice of the diagram, and obtain a set of non-rejected true effect sizes.
<p><br><b>Confidence intervals, and lower and upper confidence bounds</b>
<p>The first diagram shows the smallest possible non-rejection region. This is the diagram for computing confidence intervals. Looking at each vertical slice, the non-rejection region is centered on the point where the estimate is equal to the true value. It covers the densest region, so it can be quite compact and still contain the estimate with probability 1-α. Now looking at each horizontal slice, we see that we will always obtain non-rejection regions centered on the estimate.
<p>The second diagram shows a procedure for obtaining a lower confidence bound on the effect. Looking at each vertical slice, the non-rejection region is no longer centered -- it goes down to negative infinity. However this means the top of the non-rejection region can be moved down slightly compared to the confidence interval diagram and still contain the estimate with probability 1-α. Now looking at each horizontal slice, we see that we will obtain a non-rejection region from slightly below the estimate up to positive infinity, so this diagram is for giving lower bounds on the true effect. It gives a slightly tighter lower bound than the confidence region diagram.
<p>In a similar way we may obtain an upper confidence bound. <font size=-2>(However it would not be valid to apply both the lower bound and upper bound procedure to the same data -- in doing so we would risk rejecting the truth with probability 2α. So in this case we would need to use the confidence interval procedure.)</font>
<p><br><b>t-tests to determine the sign of an effect</b>
<p>The t-test only tests whether a true effect of zero can be rejected. Having performed this test, what can we say about other effect values?
<p>We have only tested our estimate against the boundaries at a true effect of zero, so we have only compared our estimate to the boundaries in this vertical slice. This splits the confidence interval diagram into three layers, and the confidence bound diagrams into two layers each. Looking at vertical segments within these layers, if any point is non-rejected all must be non-rejected. Filling out the confidence interval diagram and bound diagrams in this way gives us the diagrams for the two-sided and one-sided t-tests.
<p>Examining horizontal lines through these diagrams, possible outcomes are that non-rejected effects are restricted to positive values, or to negative values, or that no rejection occurs. So when the two-sided t-test rejects a true effect of zero it also tells us the sign of the effect.
<p>Similar to the lower and upper bounds, one-sided t-tests need a smaller estimate than the two-sided test to determine if the effect has a certain sign, at the cost of not testing the opposite sign.
<p>From these diagrams we can see that confidence intervals and bounds tell us more than the corresponding t-test. One small virtue of the t-test is that when it is reported a p-value for the test can be quoted, allowing the reader to set their own α and rejecting an effect of zero if p ≤ α. Confidence intervals and bounds need α to be specified before performing the procedure.
<p><br><b>TREAT</b>
<p><a href="https://academic.oup.com/bioinformatics/article-lookup/doi/10.1093/bioinformatics/btp053">The TREAT procedure</a>
(McCarthy and Smyth, 2009) represents a third kind of statistical procedure, a blend between the confidence interval and lower and upper bounds. The authors apply this to microarray and RNA-Seq data analysis, and it has implementations in their <a href="http://www.bioconductor.org/packages/release/bioc/html/limma.html">limma</a>
and
<a href="http://www.bioconductor.org/packages/release/bioc/html/edgeR.html">edgeR</a>
packages. However there is nothing stopping it from being used more generally.
<p>As described in the paper, the TREAT procedure calculates a test statistic and from this a p-value. Translating this into diagram-form: For each true effect size, a non-rejection interval centered on <i>zero</i> is found. For effect sizes close to zero, this interval is similar to our confidence interval, but further from zero it resembles an upper or lower bound.
<p>The acceptance intervals are not centered on the true effect as in confidence intervals, but are also not infinitely off-center as in the upper and lower bounds.
<p>Looking at horizontal slices through this diagram, we see that absolute effect sizes smaller than a certain amount will be rejected. This is what TREAT does: shows that the effect size is larger than a specified amount.
<p>(As TREAT is implemented in limma and edgeR, one specifies a minimum effect size and obtains a p-value. Similar to our conversion from confidence interval to t-test, this would give us a squared up H-shaped region.)
<p>Similarly to the t-test, having obtained statistical significance by TREAT, can we say anything about the sign of the effect? It seems not. There is no horizontal line which entirely rejects one or other sign. However it would be possible to fix this. This is the final diagram, my proposal for a modified TREAT, in which you learn not only that the absolute effect size is larger than some amount, but also whether it is positive or negative. Looking vertically, the non-rejection intervals are about half as far off-center. Looking horizontally, we are now able to determine the sign.
<p>The plot below shows this with t-distributed errors (df=5) and α=0.05. The diagonal lines show the boundaries for confidence intervals, lower and upper bounds, and the central line of estimate=true. You can see there is only a tiny loss of power from this modification. <font size=-2>(One would not reasonably attempt statistical analysis with a df smaller than this. For higher df the difference becomes even smaller, but does not entirely disappear.)</font>
<p><center><img src="/pfh-files/blog/01491087484/plot.svg" height=400></center>
<p><b>Conclusion</b>
<p>I hope these diagrams have given you a clearer understanding of some commonly used classical statistical procedures. They've certainly been necessary for me in order to think clearly about the TREAT procedure.
<p><br><b>References</b>
<p>McCarthy, D. J., and Smyth, G. K. (2009). Testing significance relative to a fold-change threshold is a TREAT. <i>Bioinformatics</i>, 25(6), 765-771.
Finding your ut-re-mi-fa-sol-la on a monochord, and making simple drinking straw reedpipes2017-03-16T23:29:51Z2017-03-16T23:29:51Zhttp://www.logarithmic.net/pfh/blog/01489706991
<p>I gave a class at Rowany Festival in 2013 on Guido of Arezzo's music textbook "Micrologus". It also covered a little bit on reedpipes. I don't seem to have posted the handout online. Rectifying that now:
<p>* <a href="http://www.logarithmic.net/pfh-files/blog/01489706991/monochords_and_reedpipes.pdf">Handout (PDF)</a>
<p>* <a href="http://www.logarithmic.net/pfh-files/blog/01489706991/reedpipes2.svg">Reedpipe template (SVG)</a>Briefly,2016-12-07T00:50:28Z2016-12-07T00:50:28Zhttp://www.logarithmic.net/pfh/blog/01481071828
<p>
<center style="font-size: 400%">MAKE<br><strike>WORK</strike></center>
Shiny interactivity with grid graphics2016-10-18T05:53:19Z2016-10-18T05:53:19Zhttp://www.logarithmic.net/pfh/blog/01476769999
<p>This post may be of interest to a very small number of people. Using a fancy <a href="https://stat.ethz.ch/R-manual/R-devel/library/grid/html/grid-package.html"><span class="mono">grid</span></a>-based graphics library in R such as
<a href="https://cran.r-project.org/web/packages/lattice/index.html"><span class="mono">lattice</span></a> and want some
<a href="http://shiny.rstudio.com/articles/plot-interaction.html">Shiny brushing goodness</a>? Sadly, Shiny only supports base graphics and ggplot2.
<p>Solution: First work out the name of the relevant viewport (<span class="mono">grid.ls(viewports=T,grobs=F)</span>), then use the <a href="https://cran.r-project.org/web/packages/gridBase/index.html"><span class="mono">gridBase</span></a> package to set up a corresponding base graphics viewport.
<p><pre>
output$plot <- renderPlot({
# ... draw some grid graphics ...
seekViewport("...someviewport...")
par(new=T, plt=gridBase::gridPLT())
plot(1, type="n", axes=F,
xlab="", ylab="",
xlim=c(0,1), ylim=c(0,1),
xaxs="i", yaxs="i")
})
</pre>
<p>If the viewport isn't named, tell the developer of the package to name their viewports.
Crash course in R a la 2016 with a biological flavour2016-08-07T02:17:31Z2016-08-07T02:17:31Zhttp://www.logarithmic.net/pfh/blog/01470536251
<p><ul><li><a href="https://monashbioinformaticsplatform.github.io/r-more/">Course notes</a></li></ul>
<p>Monash Bioinformatics Platform gave this day long course on Friday 5 July, 2016 (following an <a href="https://monashbioinformaticsplatform.github.io/r-intro/">introductory course</a> on Thursday). The course style is inspired by Software Carpentry, but currently sacrificing some interaction to get where I wanted to go, which is someone able to do useful bioinformatics work at a grad-student level.
<p>We definitely need to work on getting more interaction going, everyone was too quiet. This may have just been an effect of packing in so much material. However, plenty of solutions to challenges via the chat system, fairly sure most people were more or less following. The anonmymity of chat introduces an interesting dynamic -- students can propose solutions without worrying about defending their ego or propriety, or worrying about stereotype threat.Vectors are enough2016-07-17T01:45:24Z2016-07-17T01:45:24Zhttp://www.logarithmic.net/pfh/blog/01468719924
Numpy goes to some lengths to support multi-dimensional arrays, "tensors". This isn't necessary, vectors are enough.
<p>R has a different solution. A data frame (similar to a table in a database) with column vectors for the dimensions and for the value stores the same information, and a data frame can store many other things besides. Sparse tensors are easy, for example.
<p>As Judea Pearl might say, statistics is not about causality. A tensor implies causality, the dimensions of the tensor <i>cause</i> the values it holds. R, being a statistical language, rightly does not require causality to be specified and gains flexibility. Or we might say that the data frame approach does not require the primary key columns to be made explicit. So for example a bijective map need not privilege one direction of mapping (as a Python dict would).
<p>A data frame dense in some of its columns (eg from a full factorial experiment) can be stored more compactly as a tensor, but this is a detail of implementation that could be hidden (eg in a new flavour of <a href="https://cran.r-project.org/web/packages/dplyr/vignettes/databases.html">tibble</a>). Similarly, one might want various indexes to be maintained alongside a data frame.