Appendix: Cross-cultural predictors of rape

d = na.omit(sccsfac[qw(v166, v167, v667)])
d = data.frame(
    # Recode v166 and v167 so higher values mean more premarital
    # sex.
    prem.male = 4 - as.integer(d$v166),
    prem.female = 4 - as.integer(d$v167),
    rape = d$v667 == "Present")
c("Complete cases" = nrow(d))

It turns out that prem.male and prem.female are generally equal, and neither varies much:

with(d, table(prem.male, prem.female))

Suppose we then just use prem.female.

with(d, table(prem.female, rape))

It looks like rape is about 50% true across with board with an oddball at prem.female == 2. Let's try logistic regression just for laughs.

model = glm(rape ~ prem.female, data = d,
    family = binomial(link = "logit"))

Not surprisingly, this model is completely useless. The probability it predicts of rape being true is greater than 1/2 for every case:

round(d = 3, min(predict(model, type = "response")))

Conclusion: across societies, premarital-sex frequency is not predictive of rape.

d = na.omit(sccsfac[c(virtue, "v667")])
c("Complete cases" = nrow(d))

This number is similar to that of the previous analysis, but unfortunately, the intersection of cases is substantially smaller:

nrow(na.omit(sccsfac[c(virtue, "v166", "v167", "v667")]))

So we'll do an analysis with just the virtue variables and the DV (v667).

The virtue variables have many high correlations:

corm = cor(d[virtue], method = "spearman")
corm[!upper.tri(corm)] = NA
round(d = 2, corm)

Let's try summarizing them with a principal component or two. We'll do PCA based on the Spearman correlation matrix, since these variables are all on rating scales.

dvr = d[virtue]
for (i in seq_along(dvr))
    dvr[,i] = rank(dvr[,i])
round(d = 2, eigen(cor(dvr))$values / ncol(dvr))

The first PC accounts for the bulk of the variance (71%).

virt.pca = prcomp(dvr, center = T, scale = T)
round(d = 2, virt.pca$rotation[,1:3])

Here are loadings from the first three PCs. On the first, the loadings are such that positive values of the PC indicate more sexual restraint. The second is more opaque (its pattern is consistent with the male–female distinction, but not very strongly) and the third doesn't account for much variance (it distinguishes between the "sexual expression" and "sexual nonrestraint" variables, which the authors say should be added together anyway). So let's just use the first as a general measure of sexual restraint.

dm = data.frame(
    restraint = virt.pca$x[,"PC1"],
    rape = d$v667 == "Present")
model = glm(rape ~ restraint, data = dm,
    family = binomial(link = "logit"))

This time, at least, we have predicted probabilities on both sides of .5:

table(predict(model, type = "response") >= .5)

So how well can this model fit its own training data?

round(d = 2,
    mean((predict(model, type = "response") >= .5) == model$y))

Amazingly, the model has achieved a training-set accuracy of less than 50%.

Okay, that was a bust; how about adding PC2 for an additional 15% of the variance in the original IVs?

dm = data.frame(
    restraint.pc1 = virt.pca$x[,"PC1"],
    restraint.pc2 = virt.pca$x[,"PC2"],
    rape = d$v667 == "Present")
model = glm(rape ~ restraint.pc1 * restraint.pc2, data = dm,
    family = binomial(link = "logit"))

round(d = 2,
    mean((predict(model, type = "response") > .5) == model$y))

That's not much bigger than the base rate in the training set:

round(d = 2,
    max(mean(model$y), mean(!model$y)))

I'm not going to bother checking test error.

Conclusion: across societies, attitudes about sexual expression in children are not predictive of rape.

Look at the cross-classification table:

table(na.omit(sccsfac[qw(v626, v667)]))

This table suggests that in almost all cultures with less rape, there is no "Belief That Women Are Generally Inferior to Men". In cultures with more rape, on the other hand, such a belief is reasonably likely, although perhaps still less likely than a 50% chance.

Here is the table with probabilities of rape as a function of belief:

round(d = 2, prop.table(margin = 1,
    table(na.omit(sccsfac[qw(v626, v667)]))))

Without the belief, higher rape is 44% likely, but with it, it's 91% likely. Interesting.

c(accuracy = round(d = 2, with(na.omit(sccsfac[qw(v626, v667)]),
    mean((v626 == "Yes") == (v667 == "Present")))))

The accuracy rate is still not very impressive, since in the no-such-belief case, the rape rate is around 50%.

I won't do a separate (i.e., cross-validatory) estimate of predictive accuracy because I don't think of this model as having any free parameters. (If this IV is to be nontrivially dichotomously used to predict this DV, No such belief must map to Absent and Yes to Present. Vice versa wouldn't make sense.)

Conclusion: among cultures with a belief that women are generally inferior to men, rape is likely to be reasonably common or institutionalized. Among cultures without such a belief, rape rates vary more widely.

Here's the contingency table for female power and rape:

table(na.omit(sccsfac[qw(v663, v667)]))

An inconsistent picture. However, the modal level of v663 also seems to be one in which the rape rate is distant from 50%. So let's see how well we can do with logistic regression.

dm = with(na.omit(sccsfac[qw(v663, v667)]), data.frame(
    female.power = as.integer(v663),
    rape = v667 == "Present"))
model = glm(rape ~ female.power, data = dm,
    family = binomial(link = "logit"))

round(d = 2,
    mean((predict(model, type = "response") > .5) == model$y))

Not great.

Suppose we add ideology of male superiority. It has only three levels, so let's treat it as nominal.

dm = with(na.omit(sccsfac[qw(v663, v1767, v667)]), data.frame(
    female.power = as.integer(v663),
    male.sup = v1767,
    rape = v667 == "Present"))
model = glm(rape ~ female.power * male.sup, data = dm,
    family = binomial(link = "logit"))

round(d = 2,
    mean((predict(model, type = "response") > .5) == model$y))

Better. But all the predictive ability is coming from v1767, actually:

model = glm(rape ~ male.sup, data = dm,
    family = binomial(link = "logit"))
round(d = 2,
    mean((predict(model, type = "response") > .5) == model$y))

So let's use just that.

d = with(na.omit(sccsfac[qw(v1767, v667)]), data.frame(
    male.sup = ordered(levels = qw(None, Weak, Strong),
        ifelse(v1767 == "strongly articulated ideology of male superiority", "Strong",
        ifelse(v1767 == "weakly articulated ideology of male superiority", "Weak",
        "None"))),
    rape = v667 == "Present"))
table(d)

Here it is with probabilities:

round(d = 2, prop.table(table(d), 1))

You can see that ideologies of male superiority, weak or strong, predict rape at 70% probability, whereas societies without such an ideology are only 41% likely to have lots of rape.

round(d = 2, mean(with(d,
    (male.sup %in% qw(Weak, Strong)) == rape)))

Let's estimate predictive accuracy. The only free parameter I would consider the model to have is whether Weak implies TRUE or FALSE. None will be fixed at FALSE and Strong at TRUE.

cv.folds = 10
set.seed(20)
d$fold = sample(rep(factor(1 : cv.folds), len = nrow(d)))
results = c(recursive = T, lapply(1 : cv.folds, function(fold)
   {test = d[d$fold == fold,]
    train = d[d$fold != fold,]
    weak.rate = mean(subset(train, male.sup == "Weak")$rape)
    weak.pred =
        if (weak.rate < .5) F
        else if (weak.rate > .5) T
        else sample(c(T, F), 1)
    with(test, rape ==
        ifelse(male.sup == "None", F,
        ifelse(male.sup == "Strong", T,
        weak.pred)))}))
round(d = 2, mean(results))

The cross-validatory folds always assigned Weak to TRUE, so this is the same number as the training-set accuracy.

Conclusion: across cultures, an ideology of male superiority predicts more rape.