https://twitter.com/TimoRoettger/status/1262453590502912000

clusters.R

library(tidyverse)

# Simulate data under 1-cluster solution:
# Two Gaussians, with different means, same SD.
m1 = 10
m2 = 20
sd = 10
n1 = 50
n2 = 50
x1 = rnorm(n1, m1, sd)
x2 = rnorm(n2, m2, sd)
data = data.frame(
    group = rep(c(1, 2), c(n1, n2)),
    x = c(x1, x2))
ggplot(data, aes(x, fill=factor(group), group=group)) +
    geom_histogram(position='identity', alpha=.5, binwidth=10)

# Means per condition
data %>% group_by(group) %>% summarise_at(vars(x), funs(mean=mean, sd=sd))
## # A tibble: 2 x 3
##   group  mean    sd
##   <dbl> <dbl> <dbl>
## 1     1  9.58  11.8
## 2     2 19.4   11.6

# Fit 2-cluster solution
k = 2
km = kmeans(x, k)

if(km$centers[1] > km$centers[2]){
    data$com = ifelse(km$cluster==1, 1, 0)  # Cluster 1 is COM
} else {
    data$com = ifelse(km$cluster==2, 1, 0) # Cluster 2 is COM
}



# Proportion in COM cluster
data %>% group_by(group) %>% summarise_at(vars(com), funs(mean=mean, sd=sd))
## # A tibble: 2 x 3
##   group  mean    sd
##   <dbl> <dbl> <dbl>
## 1     1  0.38 0.490
## 2     2  0.68 0.471
glm(com ~ group, data=data, family=binomial) %>% summary()
##
## Call:
## glm(formula = com ~ group, family = binomial, data = data)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -1.5096  -0.9778   0.8782   0.8782   1.3911
##
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -1.7329     0.6569  -2.638  0.00834 **
## group         1.2433     0.4205   2.957  0.00311 **
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
##     Null deviance: 138.27  on 99  degrees of freedom
## Residual deviance: 129.09  on 98  degrees of freedom
## AIC: 133.09
##
## Number of Fisher Scoring iterations: 4