Reanalysis of Reynolds et al. (2020) 10.1073/pnas.2014403117

Reynolds2020.R

# Reanalysis of Reynolds et al. (2020) 10.1073/pnas.2014403117
# By Nick Brown, November 2020
# Licence: CC-0
# Dataset file available for download from https://www.pnas.org/content/suppl/2020/11/11/2014403117.DCSupplemental

df <- read.csv("pnas.2014403117.sd01.csv", stringsAsFactors=FALSE)

# Remove cases with missing data on the two outcome measures
df.crp <- df[!is.na(df$elevatedCRP),]
df.ht <- df[!is.na(df$htn),]

cat("Inflammation (CRP) analysis for patrilineal groups\n")
df.crp.pat <- df.crp[df.crp$Matriliny == 0,]
df.crp.pat.female <- df.crp.pat[df.crp.pat$Male == 0,]
df.crp.pat.male <- df.crp.pat[df.crp.pat$Male == 1,]

cont.crp.pat <- matrix(c(
  sum(df.crp.pat.female$elevatedCRP == 0)
, sum(df.crp.pat.female$elevatedCRP == 1)
, sum(df.crp.pat.male$elevatedCRP == 0)
, sum(df.crp.pat.male$elevatedCRP == 1)
), ncol=2)

ft <- fisher.test(cont.crp.pat)

female.cases <- cont.crp.pat[2]
female.total <- female.cases + cont.crp.pat[1]
female.perc <- sprintf("%.2f", (100 * female.cases / female.total))
cat("Females: total=", female.total, " cases=", female.cases, " (", female.perc, "%)", "\n", sep="")
male.cases <- cont.crp.pat[4]
male.total <- male.cases + cont.crp.pat[3]
male.perc <- sprintf("%.2f", (100 * male.cases / male.total))
cat("Males: total=", male.total, " cases=", male.cases, " (", male.perc, "%)", "\n", sep="")
cat("p-value for difference by Fisher's exact test=", sprintf("%.3f", ft$p.value), "\n", sep="")
cat("\n")

cat("Hypertension analysis for patrilineal groups\n")
df.ht.pat <- df.ht[df.ht$Matriliny == 0,]
df.ht.pat.female <- df.ht.pat[df.ht.pat$Male == 0,]
df.ht.pat.male <- df.ht.pat[df.ht.pat$Male == 1,]

cont.ht.pat <- matrix(c(
  sum(df.ht.pat.female$ht == 0)
, sum(df.ht.pat.female$ht == 1)
, sum(df.ht.pat.male$ht == 0)
, sum(df.ht.pat.male$ht == 1)
), ncol=2)

ft <- fisher.test(cont.ht.pat)

female.cases <- cont.ht.pat[2]
female.total <- female.cases + cont.ht.pat[1]
female.perc <- sprintf("%.2f", (100 * female.cases / female.total))
cat("Females: total=", female.total, " cases=", female.cases, " (", female.perc, "%)", "\n", sep="")
male.cases <- cont.ht.pat[4]
male.total <- male.cases + cont.ht.pat[3]
male.perc <- sprintf("%.2f", (100 * male.cases / male.total))
cat("Males: total=", male.total, " cases=", male.cases, " (", male.perc, "%)", "\n", sep="")
cat("p-value for difference by Fisher's exact test=", sprintf("%.3f", ft$p.value), "\n", sep="")
cat("\n")

cat("Inflammation (CRP) analysis for matrilineal groups\n")
df.crp.mat <- df.crp[df.crp$Matriliny == 1,]
df.crp.mat.female <- df.crp.mat[df.crp.mat$Male == 0,]
df.crp.mat.male <- df.crp.mat[df.crp.mat$Male == 1,]

cont.crp.mat <- matrix(c(
  sum(df.crp.mat.female$elevatedCRP == 0)
, sum(df.crp.mat.female$elevatedCRP == 1)
, sum(df.crp.mat.male$elevatedCRP == 0)
, sum(df.crp.mat.male$elevatedCRP == 1)
), ncol=2)

ft <- fisher.test(cont.crp.mat)

female.cases <- cont.crp.mat[2]
female.total <- female.cases + cont.crp.mat[1]
female.perc <- sprintf("%.2f", (100 * female.cases / female.total))
cat("Females: total=", female.total, " cases=", female.cases, " (", female.perc, "%)", "\n", sep="")
male.cases <- cont.crp.mat[4]
male.total <- male.cases + cont.crp.mat[3]
male.perc <- sprintf("%.2f", (100 * male.cases / male.total))
cat("Males: total=", male.total, " cases=", male.cases, " (", male.perc, "%)", "\n", sep="")
cat("p-value for difference by Fisher's exact test=", sprintf("%.3f", ft$p.value), "\n", sep="")
cat("\n")

cat("Hypertension analysis for matrilineal groups\n")
df.ht.mat <- df.ht[df.ht$Matriliny == 1,]
df.ht.mat.female <- df.ht.mat[df.ht.mat$Male == 0,]
df.ht.mat.male <- df.ht.mat[df.ht.mat$Male == 1,]

cont.ht.mat <- matrix(c(
  sum(df.ht.mat.female$ht == 0)
, sum(df.ht.mat.female$ht == 1)
, sum(df.ht.mat.male$ht == 0)
, sum(df.ht.mat.male$ht == 1)
), ncol=2)

ft <- fisher.test(cont.ht.mat)

female.cases <- cont.ht.mat[2]
female.total <- female.cases + cont.ht.mat[1]
female.perc <- sprintf("%.2f", (100 * female.cases / female.total))
cat("Females: total=", female.total, " cases=", female.cases, " (", female.perc, "%)", "\n", sep="")
male.cases <- cont.ht.mat[4]
male.total <- male.cases + cont.ht.mat[3]
male.perc <- sprintf("%.2f", (100 * male.cases / male.total))
cat("Males: total=", male.total, " cases=", male.cases, " (", male.perc, "%)", "\n", sep="")
cat("p-value for difference by Fisher's exact test=", sprintf("%.3f", ft$p.value), "\n", sep="")
cat("\n")

cat("Inflammation (CRP) analysis for females across groups\n")
df.crp.female <- df.crp[df.crp$Male == 0,]
df.crp.female.pat <- df.crp.female[df.crp.female$Matriliny == 0,]
df.crp.female.mat <- df.crp.female[df.crp.female$Matriliny == 1,]

cont.crp.female <- matrix(c(
  sum(df.crp.female.pat$elevatedCRP == 0)
, sum(df.crp.female.pat$elevatedCRP == 1)
, sum(df.crp.female.mat$elevatedCRP == 0)
, sum(df.crp.female.mat$elevatedCRP == 1)
), ncol=2)

ft <- fisher.test(cont.crp.female)

pat.cases <- cont.crp.female[2]
pat.total <- pat.cases + cont.crp.female[1]
pat.perc <- sprintf("%.2f", (100 * pat.cases / pat.total))
cat("Patrilineal: total=", pat.total, " cases=", pat.cases, " (", pat.perc, "%)", "\n", sep="")
mat.cases <- cont.crp.female[4]
mat.total <- mat.cases + cont.crp.female[3]
mat.perc <- sprintf("%.2f", (100 * mat.cases / mat.total))
cat("Matrilineal: total=", mat.total, " cases=", mat.cases, " (", mat.perc, "%)", "\n", sep="")
cat("p-value for difference by Fisher's exact test=", sprintf("%.3f", ft$p.value), "\n", sep="")
cat("\n")

cat("Inflammation (CRP) analysis for males across groups\n")
df.crp.male <- df.crp[df.crp$Male == 1,]
df.crp.male.pat <- df.crp.male[df.crp.male$Matriliny == 0,]
df.crp.male.mat <- df.crp.male[df.crp.male$Matriliny == 1,]

cont.crp.male <- matrix(c(
  sum(df.crp.male.pat$elevatedCRP == 0)
, sum(df.crp.male.pat$elevatedCRP == 1)
, sum(df.crp.male.mat$elevatedCRP == 0)
, sum(df.crp.male.mat$elevatedCRP == 1)
), ncol=2)

ft <- fisher.test(cont.crp.male)

pat.cases <- cont.crp.male[2]
pat.total <- pat.cases + cont.crp.male[1]
pat.perc <- sprintf("%.2f", (100 * pat.cases / pat.total))
cat("Patrilineal: total=", pat.total, " cases=", pat.cases, " (", pat.perc, "%)", "\n", sep="")
mat.cases <- cont.crp.male[4]
mat.total <- mat.cases + cont.crp.male[3]
mat.perc <- sprintf("%.2f", (100 * mat.cases / mat.total))
cat("Matrilineal: total=", mat.total, " cases=", mat.cases, " (", mat.perc, "%)", "\n", sep="")
cat("p-value for difference by Fisher's exact test=", sprintf("%.3f", ft$p.value), "\n", sep="")
cat("\n")

cat("Hypertension analysis for females across groups\n")
df.ht.female <- df.ht[df.ht$Male == 0,]
df.ht.female.pat <- df.ht.female[df.ht.female$Matriliny == 0,]
df.ht.female.mat <- df.ht.female[df.ht.female$Matriliny == 1,]

cont.ht.female <- matrix(c(
  sum(df.ht.female.pat$ht == 0)
, sum(df.ht.female.pat$ht == 1)
, sum(df.ht.female.mat$ht == 0)
, sum(df.ht.female.mat$ht == 1)
), ncol=2)

ft <- fisher.test(cont.ht.female)

pat.cases <- cont.ht.female[2]
pat.total <- pat.cases + cont.ht.female[1]
pat.perc <- sprintf("%.2f", (100 * pat.cases / pat.total))
cat("Patrilineal: total=", pat.total, " cases=", pat.cases, " (", pat.perc, "%)", "\n", sep="")
mat.cases <- cont.ht.female[4]
mat.total <- mat.cases + cont.ht.female[3]
mat.perc <- sprintf("%.2f", (100 * mat.cases / mat.total))
cat("Matrilineal: total=", mat.total, " cases=", mat.cases, " (", mat.perc, "%)", "\n", sep="")
cat("p-value for difference by Fisher's exact test=", sprintf("%.3f", ft$p.value), "\n", sep="")
cat("\n")

cat("Hypertension analysis for males across groups\n")
df.ht.male <- df.ht[df.ht$Male == 1,]
df.ht.male.pat <- df.ht.male[df.ht.male$Matriliny == 0,]
df.ht.male.mat <- df.ht.male[df.ht.male$Matriliny == 1,]

cont.ht.male <- matrix(c(
  sum(df.ht.male.pat$ht == 0)
, sum(df.ht.male.pat$ht == 1)
, sum(df.ht.male.mat$ht == 0)
, sum(df.ht.male.mat$ht == 1)
), ncol=2)

ft <- fisher.test(cont.ht.male)

pat.cases <- cont.ht.male[2]
pat.total <- pat.cases + cont.ht.male[1]
pat.perc <- sprintf("%.2f", (100 * pat.cases / pat.total))
cat("Patrilineal: total=", pat.total, " cases=", pat.cases, " (", pat.perc, "%)", "\n", sep="")
mat.cases <- cont.ht.male[4]
mat.total <- mat.cases + cont.ht.male[3]
mat.perc <- sprintf("%.2f", (100 * mat.cases / mat.total))
cat("Matrilineal: total=", mat.total, " cases=", mat.cases, " (", mat.perc, "%)", "\n", sep="")
cat("p-value for difference by Fisher's exact test=", sprintf("%.3f", ft$p.value), "\n", sep="")
cat("\n")