Recent content on Biased and Inefficient
Derek Lowe, writing on models for protein-ligand interaction Ideally, you train your model on a big pile o’ protein data, while holding out a good piece of it to keep it out of the training set. Then once your software has machine-learned its way to happiness, you see how it performs on those cases that it’s never seen before. But here’s the problem: it is very hard to set this up in a way so that…
The almost-sure representation theorem in probability says that if \(X_n\stackrel{d}{\to}X\) we can find a sequence \(\tilde X_n, \tilde X\) with the same distributions as \(X_n,X\), possibly on a different probability space, so that \(\tilde X_n\stackrel{a.s.}{\to}\tilde X\) (where the a.s. is with respect to the distribution of \(X\)). General versions of this theorem are a bit tricky to prove,…
As I have said before, the distinction between multistage and multiphase surveys is mostly important as a shibboleth for survey statisticians. The two=phase quantities \[\pi_i^*=P(i\text{ sampled}|\text{phase 1})P(i\in\text{phase 1})\] are not the same as the marginal sampling probabilities \[\pi_i=P(i\text{ sampled})\] but one uses them in exactly the same formulas, so the difference doesn’t…
In one of the social media discussions about causal inference the suggestion was made that predictive models are all you need: a good predictive model gives you all the conditional distributions you could want, and you don’t need any special causal inference stuff. I think there’s something to this point of view, but there are a few limitations. The first is that causal inference theory (eg causal…
Version 4.5 of survey is on CRAN now. There are a lot of little changes and a few new features. Thanks to Stas Kolenikov we have Bell-McCaffrey standard errors for svyglm. The standard svyglm standard errors are based on sums of squares of PSU-level residuals. As Bell and McCaffrey point out These sums of squares tend to be too small for two reasons: residuals are generally smaller than true…
This post is partly because I’m about to start teaching generalised linear models and partly to avoid doomscrolling Do predictive models need to be causal? At first glance the response is something like “Of course not; are you high?”. If you have a model \[\log P(Y=1)= \alpha +\beta_X X+\beta_Z Z\] then whether the model is usefully predictive has almost nothing to do with whether one of…
Attention Conservation Notice: I stole the title from Cosma Shalizi, who is a better writer than me. He meant the central limit theorem; I mean the Gauss-Markov theorem. Also, this is about multiparameter tests, which are intrinsically not very interesting In survey statistics, there are two classes of tests that are widely used for contingency tables and rather less widely used for regression…
Riffing off Andrew Gelman’s (correct) advice to just simulate, in both teaching and research, suppose I have a new robust location estimator. For a data set \(X_1,\dots, X_n\), I want to find the densest half of the distribution and take some location summary of it.1 Let’s say I want the mean. That is, I want to define a 50% subset of the data by function(x) { s<-sort(x) n<-length(s) halfn<- n%/%2…
Yes This question comes up from time to time on social media or StackExchange or email, often from reasonable people, so extra emphasis might be useful. There are two parts to the answer: Yes and If you think about it, what else could it be using? Write \(\hat\beta\) for the svyglm estimator. Theory says the estimator solves the weighted score equations \[U(\beta)=\sum_{i=1}^N \frac{R_i}{\pi_i}…
Suppose we have a prediction problem. We want to predict whether \(Y=1\) or \(Y=0\), but for nearly all the examples in our training data it turns out that \(Y=0\). Many predictive techniques, faced with data like this, will degenerate to predicting \(Y=0\) everywhere; even the more successful techniques will predict \(Y=0\) for most inputs. The negative predictive value \(P(Y=0|\hat Y=0)\) will…
This is approximately my Belz Lecture for the Statistical Society of Australia (Vic/Tas), from a couple of weeks ago. ‘Approximately’, because it’s from my notes. Kia ora tātou Thanks for the introduction and thanks to the Statistical Society of Australia for inviting me . It’s traditional in a public lecture like this to have some sort of polemic. I did a polemic on Data…
There’s increasing interest in (or, at least, coverage of) the idea that surveying large language models could substitute for surveying people. This has obvious advantages: to repurpose a Heinlein quote “it’s cheap, clean, convenient, and free of any possibility of wrongdoing–and you don’t have to go home in the cold”. The problem is I can’t see how it could possibly work. There is a…
There was a question on StackOverflow1 about continuous data with extra zeroes and comparing a two-part model to a continuous-only model using AIC. This doesn’t work, but the reason is interesting and some other things that don’t work are less obvious. For a long time I vaguely wondered about what was going on if you compare a Normal model and a binomial model using AIC. Mechanically, nothing…
Suppose you have a single transferable vote election to pick some number \(d\) of people, and that you belong to a large voting bloc of sheeplike followers or of cooperative forward-thinkers (delete as necessary). How should you all vote? I think the following is true: If there is a set of \(k\) or more candidates such that everyone in your block ranks all those candidates above any candidates not…
The phrase asymptotically unbiased has at least two uses in statistics. Some people use it to mean that sequence of estimators \(\hat\theta_n\) of a parameter \(\theta\) satisfies \[\lim_{n\to\infty}E[\hat \theta_n]=\theta.\] I don’t think this is a very interesting property: it’s weaker than consistency plus existence of the first moment, and if you’re going to do asymptotics you should at least…
If you have two regression models \[E[Y|X]=\beta_0+\beta_1X\] and \[E[Y|X,Z]=\gamma_0+\gamma_1X+\gamma_2Z\] then typically \(\gamma_1\neq\beta_1\), because they are different things1 A common name for this phenomenon is omitted-variable bias. That’s an unfortunate name, because it implies a direction in a situation that’s completely symmetric. Yes, \(\hat\beta_1\) is biased for \(\gamma_1\), but…
The 2025 Ihaka Lectures are coming, in September this year. The theme is “Patterns in Space”^[yes: Patterns in S P A A A C E] and we have two speakers, both talking (in part) about how modern computing environments and open-source software improve the research and communication of spatial modellers. The lectures will be on the University of Auckland in the Science building (Building…
Attention Conservation Notice: someone asked about this on social media, and I probably should have it written down for the advanced regression course, but you probably don’t need to know it. Also, if you’re an econometrician then I hate your notation and so you will probably hate mine In two-stage least squares you have an outcome \(Y\) and want to fit a causal linear regression…
People often ask what the difference is between strata and clusters in sampling. Here’s my current soundbite version Strata are a partition of the population where we hope different strata are different. We sample from all the strata in the population, so the generalisation from the sample to the population is within strata and is more accurate because of stratified sampling. Clusters are a…
The file iris_data.csv contains some data on iris photographs from unsplash.org. For each of 28 photos we have the mean luminance, chroma, and hue and the standard deviation of luminance, chroma, and hue, based on converting from RGB to polar Lab coordinates with the colorspace package. I haven’t included the actual photo files because they’re a bit big, but here are a couple of the photos at…
As everyone knows, if you want to repeat a sequence of random numbers in statistical software1 you can reset your software’s random number stream to the same position and everything will deterministically repeat2. What’s slightly more complicated is if you want ‘the same’ random numbers from a different set of parameters. Suppose you have some simulations, which might be estimating the operating…
Again with the ordinal comparisons, yes. The ‘scale of measurement’ paradigm for variables says that ordinal data are determined only up to monotone transformation, just as interval data are determined up to translation and nominal data up to reordering. This is all good for a single observation, and as I’ve argued before a problem only arises when you want to consider a distribution over multiple…
The piranha problem is both a metaphor and a set of theorems coming out of Andrew Gelman’s research group. The metaphor is of large intervention effects as piranhas that can’t be kept together in the same tank since they’d eat each other. The theorems show that a large set of intervention effects that must add up to a large total of explained variability unless they are highly correlated with each…
Quartets of datasets are trendy now, so here’s one for collinearity diagnostics. These four data sets have two candidate predictors \(x\) and \(z\) and outcome \(y\), and have very strong correlation (over 0.95) between \(x\) and \(z\). I will use the car package for variance inflation factors set.seed(2024-11-21) library(car) ## Loading required package: carData Redundancy In the first example,…
The 1986 Joint Statistical Meetings had a Data Expo where participants were given a dataset to show off visualisation approaches. The “pollen” “data” was a \(3848\times 5\) matrix purporting to be five measurements of 3848 pollen grains. In fact, the data were artificial and consisted of an Easter egg surrounded by an ellipsoidal shell. With the displays of forty years ago it was quite hard to…
The distribution of a \(n\times n\) quadratic form in Gaussian variables, \(X^TAX\) is a linear combination of \(\chi_1^2\) variables: \[Q\sim\sum_{i=1}^n\lambda_iZ^2_i\] where \(\lambda_i\) are the eigenvalues of \(A\Xi\) and \(\Xi\) is the covariance matrix of \(X\). This is not in the tables at the back of the book. For reasonable values of \(n\) and reasonable tail probabilities, the…
The National Land Transport Plan has just been released. Greater Auckland have discussion here including a chart comparing regions. I’m going to see how hexmaps created with the DHBins package compare for inter-region comparisons. It’s a little messy because the NLTP doesn’t quite report by regions: it combines some regions First, the data, in millions of dollars projected expenditure over three…
They’re back, again. For the 2024 Ihaka Lectures we have 12 September: Yihui, who needs no introduction, on Reimagining literate programming and automated report generation 19 September: Farah Hancock, a data journalist. Initially trained as a designer, Farah Hancock worked for many years in the advertising industry in New Zealand and abroad, primarily working on digital campaigns. Experiencing a…
Attention Conservation Notice: Do you really care about multiparameter tests? Why, though? In classical statistics there are three types of test for parametric models: likelihood ratio tests, Wald tests, and score tests. These are asymptotically equivalent for local alternatives and have (asymptotically) a \(\chi^2_q\) distribution when testing a \(q\)-dimensional restriction on the parameters.…
One of the projects from the R Spring and R Dev Day has been to add unequal sampling without replacement with specified marginal sampling probabilities to R’s sample(). What R currently has is sequential sampling: if you give it a set of priorities \(w\) it will sample an element with probability proportional to \(w\) from the population, remove it from the population, then sample with probability…
The term “t-test” is a bit of a troll here, since I don’t mean either a test or the Normal-based inference developed by Gossett. I’m interested in two-sample comparisons of means, effectively non-parametric in moderate to large samples. A frequentist in an elementary course would do this by saying \(\bar X\) and \(\bar Y\) are each (roughly) Normal by the Central Limit Theorem, so \(\bar X-\bar…
When two-phase study designs started being used in epidemiology and biostatistics there was a period of conflict. Survey statisticians insisted on the term “two-phase” and biostatisticians (following survey textbooks in some cases) wanted to call these “two-stage” designs. Like the correct pronounciation of ‘Scheveningen’1, the terminology identified communities. In a \(K\)-stage survey design we…
In sandwich variance estimators, the middle of the sandwich1 is the variance of the estimating function. If we have independent observations and estimate by solving \[\sum_i U_i(\theta)=0\] we want \(\mathrm{var}[\sum_i U_i(\theta)]\), which we estimate by \(\sum_i U_i(\hat\theta)U_i(\hat\theta)^T\). There are two issues here you might miss. First, some people (younger me, for example) worry that…
The survey package returns many results as svystat objects, which are numeric vectors with variance matrix as an attribute (and other optional attributes). Because they’re not made of magic, if you transform the point estimate the variance matrix doesn’t transform and is no longer appropriate. But what if they were made of magic? We have svycontrast to do delta-method transformations and we have…
The S3 method system only lets you dispatch methods on one argument of the generic. Most people use the first argument, and it’s not unheard of for people to claim that only the first argument is allowed. Actually, other arguments can be used! What’s more, if you write functions using the old-school formula/data structure, there’s a genuine reason to dispatch on the second argument. Let’s look at…
Background The current development of the survey package now has an experimental implementation of cross-validation using replicate-weight decompositions of the data. This is experimental. It is liable to change, and may contain nuts. The basic idea, as studied by Amaia Iparraguirre is to decompose survey data in ways that respect the structure of the sampling1. Complex survey data typically have…
In dual-frame sampling you take two samples from overlapping sampling frames and you need to downweight people who could have been chosen in either frame so the overlap of the two frames isn’t counted twice. Suppose you have some constant value \(\theta\) to do the downweighting, so that people in the overlap who were sampled from frame \(A\) get their weight multiplied by \(\theta\) and people in…
I’ve mentioned before that mathematician Tim Gowers had run a ‘polymath’ (massively collaborative maths research) project on non-transitive dice. There’s an arXiv preprint. There’s also a detailed write-up in Quanta, which is a magazine devoted to popular explanations of maths. As I’ve said before, this is statistically interesting (as well as being just interesting) because any instance of…
[Updated: ran it with the right version of the code] I’m writing code for multiple-frame surveys in the survey package now, and it’s at the stage where the basic stuff works (Revision 337 from r-forge) though there’s quite a bit more to implement. The canonical references are the papers by Lohr and Rao. This post is just me thinking about it. If I had actual artistic talent or an ethically-trained…
When I wrote about weights I mentioned that there was in some senses a fourth type of weights after sampling weights, precision weights, and frequency weights. The idea is that sometimes you have weights that you want to apply to an estimating function, but that they don’t have the same ontological commitments that any of the the three sets of weights come with. I’m working on dual-frame (and…
One problem in teaching statistics and communicating statistics and so on is “assumptions”. In fact, there’s at least two problems: Necessary vs sufficient The first problem is in maths communication. In maths you write down some assumptions and show they imply a conclusion. It’s usual in statistics for the assumptions to be sufficient for the conclusion, but pretty unusual for them to be…
Two interesting examples from my e-bike: The first is common in e-bikes. The charge indicator is a set of five little rectangles inside a battery outline, which makes sense. It’s very non-linear, though. The first little rectangle is almost half the battery charge. When I mentioned this on Twitter some years ago the response was that lithium batteries are non-linear and there’s nothing that can be…
‘Symbolically nested’ is a phrase I invented to distinguish two different types of nested model when writing my book about survey analysis. There has been at least one question on Stack Overflow about the phrase, so I think it’s worth explaining in a bit more detail. Often, in math-stat discussions of nested models, you see the smaller model written with predictors \(X\) and coefficients…
With the long-awaited demise of stringsAsFactors=TRUE it’s now easier to use text strings in R. It’s good that strings don’t automatically get turned into factors at read time, but the price is that strings don’t automatically get turned into factors at read time: if you have variables that need to be factors, you have to turn them into factors yourself. Factors are still a important data type in…
If we have a domain or subpopulation \({\cal D}\) and want to estimate the mean of a variable \(Y\) in that domain, the usual survey estimator is \[\hat \mu_{\cal D}=\frac{\sum_{R_i=1} w_i Y_i I(i\in {\cal D})}{\sum_{R_i=1} w_i I(i\in {\cal D})}.\] That is, it’s the estimated population total in the domain divided by the estimated population count in the domain. We’ll call this a direct estimator;…
Version 4.4-1 of the survey package for R is percolating through CRAN. There are some important additions, visible and invisible The main invisible addition is from Ben Schneider, who has written a set of C++ routines that do the multistage stratified variance calculations previously done by svyrecvar. The compiled versions are the default; use options(survey.use_rcpp=FALSE) to disable them. The…
The DOOR outcome strategy – “Desirability Of Outcome Ranking” – is a relatively new approach to composite outcomes in clinical trials. Rather than collapsing multiple outcomes – death, heart attack, new-onset angina, bad hair day – into a single binary ‘bad thing’, the idea is to rank the trial participants by how bad their outcome is. DOOR is obviously attractive: these bad events are not all…
Suppose you find, by analysis of crime data, that a house which has been burgled once is more likely to be burgled again in the following few months. This could happen because the house just has high burglary risk, due to the neighbourhood it’s in, the availability of easy escape routes, indicators of wealth, and so on. It could also happen because burglars know to come back a few months later…
Attention Conservation Notice: This is probably well known in some circles Suppose you have a (parametric, Normal) linear mixed model \[Y=X\beta+Zb+\epsilon\] where \(\epsilon\) are iid \(N(0,\sigma^2)\) and \(b\) are \(N(0, \sigma^2V(\theta))\). Write \(\Xi\) for the marginal covariance matrix of \(Y\): \[\Xi = \mathrm{cov}[Y]=\sigma^2(I+Z^TVZ)\] The loglikelihood can be written…
Attention Conservation Notice: this is about multiparameter hypothesis tests, which are intrinsically not very interesting In regression models (and contingency tables) for survey data, there are two classes of tests based on a division that’s more or less orthogonal to the score/Wald/LRT division. Consider score tests, and for notational simplicity pretend that we’re interested a test of the…