# Power of sibling estimate of for rare variants

Author

Gibran Hemani

Published

November 10, 2023

## Background

Ascertaining families who have rare disease can improve power of detection because rare variant is segregating within the family. What about rare variants influencing trait where families are not ascertained?

``library(dplyr)``
``````
Attaching package: 'dplyr'``````
``````The following objects are masked from 'package:stats':

filter, lag``````
``````The following objects are masked from 'package:base':

intersect, setdiff, setequal, union``````
``````make_families <- function(af, nfam) {
nsnp <- length(af)
mums <- matrix(0, nfam, nsnp)
sibs1 <- matrix(0, nfam, nsnp)
sibs2 <- matrix(0, nfam, nsnp)
ibd <- matrix(0, nfam, nsnp)
ibs <- matrix(0, nfam, nsnp)
for(i in 1:nsnp)
{
dad1 <- rbinom(nfam, 1, af[i]) + 1
dad2 <- (rbinom(nfam, 1, af[i]) + 1) * -1
mum1 <- rbinom(nfam, 1, af[i]) + 1
mum2 <- (rbinom(nfam, 1, af[i]) + 1) * -1

dadindex <- sample(c(TRUE, FALSE), nfam, replace=TRUE)

mumindex <- sample(c(TRUE, FALSE), nfam, replace=TRUE)
mumh <- rep(NA, nfam)
mumh[mumindex] <- mum1[mumindex]
mumh[!mumindex] <- mum2[!mumindex]

dadindex <- sample(c(TRUE, FALSE), nfam, replace=TRUE)

mumindex <- sample(c(TRUE, FALSE), nfam, replace=TRUE)
mumh <- rep(NA, nfam)
mumh[mumindex] <- mum1[mumindex]
mumh[!mumindex] <- mum2[!mumindex]

ibd[,i] <- (as.numeric(sib1[,1] == sib2[,1]) + as.numeric(sib1[,2] == sib2[,2])) / 2

sibs1[,i] <- rowSums(abs(sib1) - 1)
sibs2[,i] <- rowSums(abs(sib2) - 1)
mums[,i] <- mum1 - 1 + abs(mum2) - 1

# l[[i]] <- (sum(sib1[,1] == sib2[,1]) / nsnp + sum(sib1[,2] == sib2[,2]) / nsnp) / 2

}

# This may not be correct - getting some really large values
ibs <- scale(sibs1) * scale(sibs2)

# Just count how many alleles are in common
ibs_unw <- abs(abs(sibs1 - sibs2) - 2) / 2
}

makePhen <- function(effs, indep, vy=1, vx=rep(1, length(effs)), my=0) {
if(is.null(dim(indep))) indep <- cbind(indep)
stopifnot(ncol(indep) == length(effs))
stopifnot(length(vx) == length(effs))
cors <- effs * vx / sqrt(vx) / sqrt(vy)
stopifnot(sum(cors^2) <= 1)
cors <- c(cors, sqrt(1-sum(cors^2)))
indep <- t(t(scale(cbind(indep, rnorm(nrow(indep))))) * cors * c(vx, 1))
y <- drop(scale(rowSums(indep)) * sqrt(vy)) + my
return(y)
}

chooseEffects <- function(nsnp, totvar, sqrt=TRUE) {
eff <- rnorm(nsnp)
aeff <- abs(eff)
sc <- sum(aeff) / totvar
out <- eff / sc
if(sqrt)
{
out <- sqrt(abs(out)) * sign(out)
}
return(out)
}

make_phenotypes <- function(fam, eff_gx, eff_xy, vx, vy, mx, my) {
lapply(fam, function(g)
{
u <- rnorm(nrow(g))
x <- makePhen(c(eff_gx), cbind(g), vy=vx, my=mx)
y <- makePhen(c(eff_xy), cbind(x), vy=vy, my=my)
return(data.frame(x=x, y=y))
})
}

join_populations <- function(l) {
mums <- do.call(rbind, lapply(l, function(x) x\$mums))
sibs1 <- do.call(rbind, lapply(l, function(x) x\$sibs1))
sibs2 <- do.call(rbind, lapply(l, function(x) x\$sibs2))
ibd <- do.call(rbind, lapply(l, function(x) x\$ibd))
ibs <- do.call(rbind, lapply(l, function(x) x\$ibs))
}

sample_populations <- function(l, n) {
index <- sort(sample(1:x, n, replace=FALSE))
l\$mums <- l\$mums[index,]
l\$sibs1 <- l\$sibs1[index,]
l\$sibs2 <- l\$sibs2[index,]
l\$ibd <- l\$ibd[index,]
l\$ibs <- l\$ibs[index,]
l\$ibs_unw <- l\$ibs_unw[index,]
return(l)
}

makephen <- function(b, g, vx=1) {
pred <- b * g
e <- rnorm(length(g), 0, sqrt(1-var(pred)))
return(pred+e)
}

# Dynastic effects
dynastic_phen <- function(fam, eff_gx, eff_xy, eff_ux, eff_uy, eff_xu) {
n <- nrow(fam\$sibs1)

# parents x

# umums <- rnorm(n)
# xmums <- makePhen(c(eff_gx, eff_ux), cbind(fam\$mums, umums))

x1 <- makephen(eff_gx, fam\$sibs1)
x2 <- makephen(eff_gx, fam\$sibs2)
l <- list()
# l\$mums <- data.frame(x=xmums)
l\$sibs1 <- data.frame(x=x1)
l\$sibs2 <- data.frame(x=x2)
return(l)
}

sibreg <- function(g1, g2, y1, y2) {
d <- bind_rows(
tibble(
y=y1, fg=(g1+g2)/2, cg=g1-fg
),
tibble(
y=y2, fg=(g1+g2)/2, cg=g2-fg
)
)
m <- summary(lm(y ~ cg + fg, data=d))\$coef
tibble(bhat=m[2,1], se=m[2,2], pval=m[2,4], n=nrow(d), method="sibreg")
}

popreg <- function(y, x) {
index <- is.finite(y) & is.finite(x)
n <- sum(index)
y <- y[index]
x <- x[index]
vx <- var(x)
bhat <- cov(y, x) / vx
ahat <- mean(y) - bhat * mean(x)
# fitted <- ahat + x * bhat
# residuals <- y - fitted
# SSR <- sum((residuals - mean(residuals))^2)
# SSF <- sum((fitted - mean(fitted))^2)

rsq <- (bhat * vx)^2 / (vx * var(y))
fval <- rsq * (n-2) / (1-rsq)
tval <- sqrt(fval)
se <- abs(bhat / tval)

# Fval <- (SSF) / (SSR/(n-2))
# pval <- pf(Fval, 1, n-2, lowe=F)
p <- pf(fval, 1, n-2, lowe=F)
return(tibble(
bhat=bhat, se=se, fval=fval, pval=p, n=n, method="popreg"
))
}

sibreg2 <- function(fam, phen) {
sdiffgx <- fam\$sibs1[,1] - fam\$sibs2[,1]
sdiffx <- phen\$sibs1\$x - phen\$sibs2\$x
popreg(sdiffx, sdiffgx) %>% mutate(method="sibreg")
}

reg <- function(f, p) {
# sibreg(f\$sibs1[,1], f\$sibs2[,1], p\$sibs1\$x, p\$sibs2\$x)
bind_rows(
sibreg2(f, p),
popreg(p\$sibs1[,1], f\$sibs1[,1])
)
}

r2 <- 0.001
param <- tibble(
af=seq(0.0001, 0.01, by=0.00005),
vg=r2,
varg=2*af*(1-af),
b=sqrt(r2/varg)
)

res <- lapply(1:nrow(param), \(i) {
f <- make_families(param\$af[i], 500000)
p <- dynastic_phen(f, param\$b[i], 0, 0, 0, 0)
r <- reg(f, p)
r\$af=param\$af[i]
return(r)
}) %>% bind_rows()``````
``````library(ggplot2)
res %>% group_by(af) %>%
summarise(fvalratio=fval[1]/fval[2]) %>%
ggplot(., aes(x=af, y=fvalratio)) +
geom_point()``````

## Summary

• Power ratio does not change by allele frequency

``sessionInfo()``
``````R version 4.3.2 (2023-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

time zone: Europe/London
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] ggplot2_3.4.2 dplyr_1.1.4

loaded via a namespace (and not attached):
[1] vctrs_0.6.4       cli_3.6.1         knitr_1.45        rlang_1.1.2
[5] xfun_0.41         generics_0.1.3    jsonlite_1.8.7    labeling_0.4.2
[9] glue_1.6.2        colorspace_2.1-0  htmltools_0.5.7   scales_1.2.1
[13] fansi_1.0.5       rmarkdown_2.25    grid_4.3.2        munsell_0.5.0
[17] evaluate_0.23     tibble_3.2.1      fastmap_1.1.1     yaml_2.3.7
[21] lifecycle_1.0.4   compiler_4.3.2    htmlwidgets_1.6.3 pkgconfig_2.0.3
[25] farver_2.1.1      digest_0.6.33     R6_2.5.1          tidyselect_1.2.0
[29] utf8_1.2.4        pillar_1.9.0      magrittr_2.0.3    withr_2.5.2
[33] tools_4.3.2       gtable_0.3.3     ``````