Arnold Kling: 2021 book titles show epistemological crisis.

This may fit with a historical pattern. The barbarians sack the city, and the carriers of the dying culture repair to their basements to write.

Though some of this goes back at least to Paul Meehl.

HTT: Bob Horn. Again.

Jacobs' Heather Wishart interviews Gen. McChrystal on teams and innovation:

I think our mindset should often be that while we don’t know what we’re going to face, we’re going to develop a team that is really good and an industrial base that is fast. We are going to figure it out as it goes, and we’re going to build to need then. Now, that’s terrifying.

On military acquisition:

The mine-resistant vehicles were a classic case; they didn’t exist at all in the US inventory; we produced thousands of them. [after the war started]

Pandas is wonderful, but Haki Benita reminds us it’s often better to aggregate in the database first:

Dan Miessler reflects on the value of shared culture:

I think young men lacking a shared culture with those around them constitute a national security threat. They are empty vessels waiting to be filled with extremist ideologies around race, religion, and immigration.

(Scroll to end, look for “Declining Religion in the West…")

David Peterson suggests metascience drop its implicit unity-of-science approach:

If reformers want to make diverse scientific fields more robust, they need to demonstrate that they understand how specific fields operate before they push a set of practices about how science overall should operate.

Background

This was inspired by Paul Harrison’s (pfh’s) 2021 post, We’ve been doing k-means wrong for more than half a century. Pfh found that the K-means algorithm in his R package put too many clusters in dense areas, resulting in worse fits compared with just cutting a Ward clustering at height (k).

tl;dr

I re-implemented much of pfh’s notebook in Python, and found that Scikit-learn did just fine with k-means++ init, but reproduced the problem using naive init (random restarts). Cross-checking with flexclust, he decided the problem was a bug in the LICORS implementation of k-means++.

Upshot: use either Ward clustering or k-means++ to choose initial clusters. In Python you’re fine with Scikit-learn’s default. But curiously the Kward here ran somewhat faster.

Update Nov-2022: I just searched for the LICORS bug. LICORS hasn’t been maintained since 2013, but it’s popular in part for its implementation of kmeanspp , compared to the default (naive) kmeans in R’s stats package. However, it had a serious bug in the distance matrix computation reported by Bernd Fritzke in Nov. 2021 that likely accounts for the behavior Paul noticed. Apparently fixing that drastically improved its performance.

I’ve just created a pull-request to the LICORS package with that fix. It appears the buggy code was copied verbatim into the motifcluster​ package. I’ve added a pull request there.

Summarizing pfh

• k-means tries to minimise sum-of-squared distances of points to clusters
• The usual Lloyd’s algorithm is fast, but gets stuck near the initial guess.
• Initializing by random choice distributes clusters proportional to point density.
• But cluster centers should be more spread out, roughly as density to the power $\frac{d}{d+2}$, where $d$ is number of dimensions.
• k-means++ and Ward both seem to do this.
• In one dimension, fast exact solutions are possible. R:Ckmeans.1d.dp, Python:ckwrap, ckmeans

References

Citation for fastcluster:

Daniel Müllner, fastcluster:Fast Hierarchical, Agglomerative Clustering Routines for R and Python, Journal of Statistical Software 53 (2013), no. 9, 1–18, URL http://www.jstatsoft.org/v53/i09/

Speed depends on many things.

• Note the up-to-1000x improvements in this 2018 discussion of sklearn’s k-means by using Cython, OpenMP, and well-vectorized code. (The ticket is closed with a [MRG] note, so may have been incorporated already.)
• And this 2020 post says Facebook’s faiss library is 8x faster and 27x more accurate than sklearn, at least on larger datasets like MNIST.
• Or a 2017 post for CFXKMeans promising 78x speedup – though note sklearn may have gotten much faster since then.

Test Times & Scores

I’ll omit the code running the tests. Defined null_fit() , do_fits(), do_splice(), functions to run fits and then combine results into a dataframe.

Borrowing an example from pfh, we will generate two squares of uniform density, the first with 10K points and the second with 1K, and find $k=200$ means. Because the points have a ratio of 10:1, we expect the ideal #clusters to be split $\sqrt{10}:1$.

Python

We see the same thing using vanilla k-means (random restarts), but the default k-means++ init overcomes it.

The Data:

SciKit KMeans: Null & Full (Naive init):

SciKit KMeans: Elkan & Kmeans++:

Ward & Polish:

More data: vary $n1, k$

What we’re seeing above is that Ward is fast and nearly as good, but not better.

Let’s collect multiple samples, varying $k$ and $n1$ as we go.

Plot Scores

We will see that KWard+Polish is often competitive on score, but seldom better.

• rows = $k$
• columns = $n1$ (with $n2 = 10 \times n1$

Pfh’s example was for $n1 = 1000$ and $k = 200$.

Plot CPU Times

Remember that polish() happens after the Ward clustering, so you should really add those two columns. But in most cases it’s on the order of the KNull.

Even with the polish step, Ward is generally faster, often much faster. The first two has curious exceptions for $n1 = 100, 1000$. I’m tempted to call that setup overhead, except it’s not there for $n1 = 10$, and the charts have different orders of magnitude for the $y$ axis.

Note that the wall-time differences are less extreme, as KMeans() uses concurrent processes. (That helps the polish() step as well, but it usually has few iterations.)

Fair enough: On uniform random data, Ward is as fast as naive K-means and as good as a k-means++ init.

How does it do on data it was designed for? Let’s make some blobs and re-run.

Test with Real Clusters

OK, so that’s how it performs if we’re just quanitizing uniform random noise. What about when the data has real clusters? Saw blob generation on the faiss example post.

Preview of some of the data we’re generating:

Scores

Ward and Polish score on par with KMeans++. For some combinations of $n1$ and $k$ the other algorithms are also on par, but for some they do notably worse.

CPU Times

Ward is constant for a given $n1$, about 4-5s for $n1 = 10,000$. KMeans gets surprisingly slow as $k$ increases, taking 75s vs. Ward’s 4-5s for $k=200$. Even more surprising, Elkan is uniformly slower than full. This could be intepreted vs. compiled, but needs looking into.

The basic result holds.

Strangely, Elkan is actually slower than the old EM algorithm, despite having well-organized blobs where the triangle inequality was supposed to help.

Conclusion

On both uniform random data and blob data, KWard+Polish scores like k-means++ while running as fast as vanilla k-means (random restarts).

In uniform data, the polish step seems to be required to match k-means++. In blobs, you can pretty much stop after the initial Ward.

Surprisingly, for sklearn the EM algorithm (algorithm='full') is faster than the default Elkan.

Appendix: the Classes

We defined two classes for the tests:

• KNull - Pick $k$ data points at at random and call that your answer.
• KWard - A drop-in replacement for Scikit-Learn’s KMeans, but uses Ward hierarchical clustering, but at height $k$.

We call the fastcluster package for the actual Ward clustering, and provide a .polish() method to do a few of the usual EM iterations to polish the result.

It gives results comparable to the default k-means++ init, but (oddly) was notably faster for large $k$. This is probably just interpreted versus compiled, but needs some attention. Ward is $O(n^2)$ while k-means++ is $O(kn)$, but Ward was running in 4-5s while scikit-learn’s kmeans ++ was taking notably longer for $k≥10$. For $k=200$ it took 75s. (!!)

The classes are really a means to an end here. The post is probably most interesting as:

• identifying a bug in an R package,
• learning a bit more about what makes K-means solutions good
• finding an odd slowness in scikit-learn’s k-means.

KNull Class

KNull just choses $k$ points at random from $X$. We could write that from scratch, but it’s equivalent to calling KMeans with 1 init and 1 iteration. (Besides, I was making a subclassing mistake in KWard, and this minimal subclass let me track it down.)

class KNull(KMeans):
    """KMeans with only 1 iteration: pick $k$ points from $X$."""
    def __init__(self, 
                 n_clusters: int=3, 
                 random_state: RandLike=None):
        """Initialize w/1 init, 1 iter."""
        super().__init__(n_clusters,
                        init="random",
                        n_init=1,
                        max_iter=1,
                        random_state=random_state)  

Aside: A quick test confirms .inertia_ stores the training .score(mat). Good because it runs about 45x faster.

KWard Class

We make this a subclass of KMeans replacing the fit() method with a call to fastcluster.linkage_vector(X, method='ward') followed by cut_tree(k) to get the initial clusters, and a new polish() method that calls regular KMeans starting from the Ward clusters, and running up to 100 iterations, polishing the fit. (We could put polish() into fit() but this is clearer for testing.

Note: The _vector is a memory-saving variant. The Python fastcluster.linkage*() functions are equivalent to the R fastcluster.hclust*() functions.

class KWard(KMeans):
    """KMeans but use Ward algorithm to find the clusters.
    
    See KMeans for full docs on inherited params.
    
    Parameters
    -----------
    These work as in KMeans:
        n_clusters : int, default=8
        verbose : int, default=0
        random_state : int, RandomState instance, default=None

    THESE ARE IGNORED: 
        init, n_init, max_iter, tol
        copy_x, n_jobs, algorithm
    
    
    Attributes 
    -----------
    algorithm : "KWard"
    
    Else, populated as per KMeans:
        cluster_centers_
        labels_
        inertia_
        n_iter_ ????
    
    Notes
    -------
    Ward hierarchical clustering repeatedly joins the two most similar points. We then cut the
    resulting tree at height $k$.  We use the fast O(n^2) implementation in fastcluster.
        
    """
    def __init__(self, n_clusters: int=8, *,
                 verbose: int=0, random_state=None):
        super().__init__(n_clusters,
                         verbose = verbose, 
                         random_state = random_state)
        self.algorithm = "KWard" # TODO: Breaks _check_params()
        self.polished_ = False
        
    def fit(self, X: np.array, y=None, sample_weight=None):
        """Find K-means cluster centers using Ward clustering.
        
        Set .labels_, .cluster_centers_, .inertia_.
        Set .polished_ = False.
        
        Parameters
        ----------
        X : {array-like, sparse matrix} of shape (n_samples, n_features)
            Passed to fc.linkage_vector.
            
        y : Ignored
        
        sample_weight : Ignored
            TODO: does fc support weights?
        
        """
        # TODO - add/improve the validation/check steps
        X = self._validate_data(X, accept_sparse='csr',
                               dtype=[np.float64, np.float32],
                               order='C', copy=self.copy_x,
                               accept_large_sparse=False)
        
        # Do the clustering. Use pandas for easy add_col, groupby.
        hc = fc.linkage_vector(X, method="ward")
        dfX = pd.DataFrame(X)
        dfX['cluster'] = cut_tree(hc, n_clusters=self.n_clusters)

        # Calculate centers, labels, inertia
        _ = dfX.groupby('cluster').mean().to_numpy()
        self.cluster_centers_ = np.ascontiguousarray(_)
        self.labels_ = dfX['cluster'] 
        self.inertia_ = -self.score(X)
        
        # Return the raw Ward clustering assignment
        self.polished_ = False
        return self
    
    def polish(self, X, max_iter: int=100):
        """Use KMeans to polish the Ward centers. Modifies self!"""
        if self.polished_:
            print("Already polished. Run .fit() to reset.")
            return self
        
        # Do a few iterations
        ans = KMeans(self.n_clusters, n_init=1, max_iter=max_iter,
                     init=self.cluster_centers_)\
                .fit(X)
        
        # How far did we move?
        𝛥c = np.linalg.norm(self.cluster_centers_ - ans.cluster_centers_)    
        𝛥s = self.inertia_ - ans.inertia_
        print(f"   Centers moved by: {𝛥c:8.1f};\n"
              f"  Score improved by: {𝛥s:8.1f} (>0 good).")
        self.labels_ = ans.labels_
        self.inertia_ = ans.inertia_
        self.cluster_centers_ = ans.cluster_centers_
        self.polished_ = True
        return self

Usage: KWard(k).fit(X) or KWard(k).fit(X).polish(X).

Alan Jacobs today:

Wondering … what to read? … Ask yourself one question: Does this writer make bank when we hate one another? And if the answer is yes, don’t read that writer.

Michele Norris describes how Germany faces its Nazi past, as a model for how America could face slavery. It turns out they have a long word for it, Vergangenheitsaufarbeitung,

an abstract, polysyllabic way of saying, ‘We have to do something about the Nazis.’

V24g took awhile to get started, and it’s messy and imperfect, but it’s real. Norris lists many examples. One:

Cadets training to become police officers in Berlin take 2½ years of training that includes Holocaust history and a field trip to the Sachsenhausen concentration camp.

Her source:

“[P]olice fatally shot 11 people and injured 34 while on duty in 2018, according to statistics compiled by the German Police Academy in Münster.

… “In Minnesota alone, where Mr. Floyd was killed, police fatally shot 13 people.”

And on the 40th anniversary of the end of WW2, West German President von Weizsäcker, son of a chief Nazi diplomat, said:

We need to look truth straight in the eye. …anyone who closes his eyes to the past is blind to the present. Whoever refuses to remember the inhumanity is prone to new risk of infection.

Measuring bot misinfo on facebook - DANMASK-19

Not exactly a surprise, but Ayers et al measure a campaign. Posts in bot-infested groups were twice as likely to misrepresent the DANMASK-19 study.

Authors suggest: penalize, ban, or counter-bot. Hm.

Sigh

Re: earlier posts:

The apparent discrepancy November/December discrepancy between higher CovidTracking counts (COVID19 deaths) and lower CDC official excess death counts has essentially vanished.

The actual excess deaths amount to +6,500 per week for December. There was just a lot of lag.

I give myself some credit for considering that I might be in a bubble, that my faith in the two reporting systems might be too strong misplaced, and looking for alternate explanations.

But in the end I was too timid in their defense: I thought only about 5K of the 7K discrepancy would be lag, and that we would see a larger role of harvesting, for example.

Periodically reread SlateStarCodex' A failure, but not of prediction. There is so much gold in there.

Finishing Off the Beaten Path by Gustavo Bondoni 📚

A delightful collection of SF/Fantasy stories with a light magical-realist feel.

Currently reading: Breaking Bread with the Dead by Alan Jacobs (µblog @ayday) 📚

Bias and Noise A clear essay by Kahneman, Sibony, & Sunstein.

Executives in the insurance company said they expected about a 10 percent difference. But the typical difference we found between two underwriters was an astonishing 55 percent of their average premium

In a long post on sustained irrationality in the markets, Vitalik Buterin describes his experience wading into to 2020 election prediction markets:

I decided to make an experiment on the blockchain that I helped to create: I bought $2,000 worth of NTRUMP (tokens that pay $1 if Trump loses) on Augur. Little did I know then that my position would eventually increase to $308,249, earning me a profit of over $56,803, and that I would make all of these remaining bets, against willing counterparties, after Trump had already lost the election. What would transpire over the next two months would prove to be a fascinating case study in social psychology, expertise, arbitrage, and the limits of market efficiency, with important ramifications to anyone who is deeply interested in the possibilities of economic institution design.

There’s a skippable technical section. His take-home is that intellectual underconfidence is a big part of why these markets can stay so wrong for so long.

But nevertheless it seems to me more true than ever that, as goes the famous Yeats quote, “the best lack all conviction, while the worst are full of passionate intensity.”

Readibility or performance? In How performance became the nemesis of secure Python code, Dima Kotik argues that Python worries too much about performance.

1. Readability is its own reward.
2. Complexity is the enemy of security.
3. Many tricks don’t help in practice.

Yesterday I commented on a cousin’s post sharing a claim about 9 reported child vaccine deaths. I looked up each death in VAERS and noted two were actually gunshot wounds, 3-5 were special cases, so only 2-4 were notably concerning. I suspect this didn’t help: she quickly deleted my comment.

David Broniatowski says I shouldn’t be surprised. He argues that both debunking and censorship are counterproductive. Remember,

Russian Twitter “troll” accounts weaponized demeaning provaccine messages as frequently as vaccine refual narratives when conducting a broad campaign to promote discord in American society.

What to do instead? The hard work of opening “collaborations with public health partners”, and especially with physicians, who are generally trusted. This is of course harder. And I’m not a physician so that’s out.

Open letter from David Broniatowski:

The vaccine rollout in the USA has slowed driven, in part, by the fact that the most eager and confident citizens have now been vaccinated. The hurdle now is no longer one of vaccine supply, but rather, demand. In two new editorials, and a podcast, all in the American Journal of Public Health, I make the case that:

1. Debunking misinformation is insufficient to convince hesitant people to vaccinate. Rather, we must listen to their concerns and communicate the gist of vaccination in a manner that accords with their values.

2. Blanket removal of online content by Facebook, Twitter, and Google/YouTube may be counterproductive, driving hesitant people to seek out information on alternative platforms. On the other hand, social media platforms are excellent tools for microtargeting and can help public health agents to reach people who are the most hesitant. We can use social media, in combination, with traditional methods, to build relationships with the most hesitant people and increase their likelihood of vaccinating.

Together, these strategies can help us cross the threshold of herd immunity to end the pandemic.

Podcast is here. [<- Link may not render, but it works. -ct]

Thanks to Mike Bishop for alerting me to Jiminez' 100-tweet thread and Lancet paper on the case for COVID-19 aerosols, and the fascinating 100-year history that still shapes debate.

Because of that history, it seems admitting “airborne” or “aerosol” has been quite a sea change. Some of this is important - “droplets” are supposed to drop, while aerosols remain airborne and so circulate farther.

But some seems definitional - a large enough aerosol is hindered by masks, and a small droplet doesn’t drop.

However, point being that like measles and other respiratory viruses, “miasma” isn’t a bad concept, so contagion can travel, esp. indoors.

VAERS Caveat

Please people, if using VAERS, go check the details. @RealJoeSmalley posts stuff like “9 child deaths in nearly 4,000 vaccinations”, but it’s not his responsibility if the data is wrong, caveat emptor.

With VAERS that’s highly irresponsible - you can’t even use VAERS without reading about its limits.

I get 9 deaths in VAERS if I set the limits to “<18”. But the number of total US vaccinations for <18 isn’t 4,000 - it’s 2.2M.

Also I checked the 9 VAERS deaths for <18:

Two are concerning because little/no risk:

5. 16yo, only risk factor oral contraceptives
6. 15yo, no known risks

Two+ are concerning but seem experimental. AFAIK the vaccines are not approved for breastfeeding, and are only in clinical trial for young children. Don’t try this at home:

2. 5mo breastmilk exposure - mom vaccinated. (?!)
3. 2yo in ¿illicit? trial? Very odd report saying it was a clinical trial but the doctors would deny that, reporter is untraceable, batch info is untraceable. Odd.
4. 1yo, seizure. (Clinical trial? Else how vaccinated?)

Two were very high risk patients. (Why was this even done?):

3. 15yo with about 25 severe pre-existing / allergies
4. 17yo w/~12 severe pre-existing / allergies

Two are clearly unrelated:

1. Error - gunshot suicide found by family, but age typed as “1.08”.
2. 17yo, firearm suicide - history of mental illness

For evaluating your risk, only the two teens would seem relevant. They might not be vaccine-related, but with otherwise no known risk, it’s a very good candidate cause.

VAERS Query

I’m not able to get “saved search” to work, so here’s the non-default Query Criteria:

• Age: < 6 months; 6-11 months; 1-2 years; 3-5 years; 6-17 years
• Event Category: Death
• Serious: Yes
• Vaccine Products: COVID19 VACCINE (COVID19)
• 
  

Group By: VAERS ID

Surprisingly, fact-checking prompted Iran’s Supreme Leader to retract his mistaken claim about Iran’s economy.

It’s probably the first time in the past three decades that Ayatolla Khamenei admitted to a mistake and corrected himself.

A glimmer of hope in a social media world that treats fact-checking as evidence for the false claim. (Though before the retraction, state media photoshopped an infographic to try to bolster the claim.)

So…

Anyone know if Biden has corrected any of his (relatively rare) 4-Pinnochio claims from his first 100 days? [List.]

Or the former US President on his 3,294 unique 4-Pinnochio claims, 94 from the first 100 days? [Full dataset here].

…we are left with the problem that… Social scientists think of themselves as explorers and they will continue to sail the world’s oceans shouting “Land!” at every mirage on the horizon even if much of the Earth has already been mapped. ~Jay Greene

Wilderness & Environmental Medicine’s editor Neil Pollock frequently writes an Editor’s note about peer review, publishing, and reproducibility. I wish he were more concise, but he often raises good points, and I expect none of my colleagues read WEM, so here’s a digest. This month’s essay is Avoiding the High Cost of Peer Review Failures.
Skipping a review of the problem, and a note about unintended consequences of open access, we get to some clear caveat emptor:

The legitimacy of journals cannot be confirmed by name or impact factor scores, and often not by promises made regarding peer-review standards…. Many predatory journals have credible and even inspiring names. They can also manufacture or manipulate impact factor scores and blatantly mislead regarding peer-review practices. [including ignoring reviews]

Caveat emptor:

Mindfulness, and more than a small degree of cynicism, is necessary to critically evaluate the legitimacy of any journal.

How you will be tempted to fail:

…getting through “peer review” with no more than trivial editorial comments may seem reasonable for the person or team thinking that their words are gold.

Being invited to review may also confer an aura of legitimacy. Such events could result in additional manuscripts being submitted to the same journal.

Stop being cats:

The inherently independent nature of researchers can lead to avoidance of conversations regarding research publication. [Discuss concerns and establish institutional guidelines to avoid being trapped by predatory journals.]

For example: [breaking his sentences into bullets]

• Did a person or team publish in such a journal inadvertently or to get around research weaknesses?
• Should full (or any) credit be given for publications in journals found to be predatory?
• Should job candidates with a history of publication in predatory journals be considered?
• Should articles published in journals employing predatory practices count in tenure packages?
• What scrutiny of the effort of flagged authors is warranted?

It’s been awhile since I read Shapin, but I’m reminded of early scientific societies and the network of trust built up by personally recommending new members. At this point I can’t see submitting to a journal that isn’t already known to my colleagues and field. But replicability indices (here | here | here…) show even that is not enough - Pollock is right that your department and institution has to ask some hard questions.

Flu has nearly disappeared worldwide durin the COVID pandemic

The distancing measures to prevent COVID are super-effective against flu. (So COVID is remarkably contagious.)

Interview in Bulletin of the Atomic Scientists about his new book 📚 Restricted Data. A thoughtful, nuanced discussion. Here are just a few:

You learn that most of what they’re redacting is really boring.

Explicit information—information you can write down—by itself is rarely sufficient for these kinds of technologies. …That isn’t saying the secrets are worthless, but it is saying that they’re probably much lower value than our system believes them to be.

Once you peel back the layer of secrecy—even in the Eisenhower years—you don’t find a bunch of angry malcontented bureaucrats on the other side. You find rich discussions about what should and shouldn’t be released. You find differences of opinion, …

I was also surprised that so many aspects of the system that we’ve come to take for granted are really determined by a tiny number of people—maybe six or seven people.

HTT Bob Horn

Congratulations to Michael Gordon et al for their paper in PLoS One!

This paper combines the results of four previous replication market studies. Data & scripts are in the pooledmaRket R package.

Key points:

• Combined, it covers 103 replications in behavioral & social sciences.
• Markets were 73% accurate, surveys a bit less.
• p-values predict original findings, though not at the frequencies you’d expect.

Enough summarizing - it’s open access, go read it! 😀🔬📄

~ ~ ~

Coda

We used this knowledge in the Replication Markets project, but it took awhile to get into print, as it does.

It should be possible to get 80-90% accuracy:

• These were one-off markets - no feedback and no learning!
• A simple p-value model does nearly as well, with different predictions.
• Simple NLP models on the PDF of the paper do nearly as well, with different predictions.

Replication Markets probably did worse ☹️, but another team may have done better. TBD.