[R] Pooled cohort equation (10-year absolute risk)

## Pooled cohort equation: implementing pooled cohort equation
## AR10 (ten-year absolute risk) = 1 - S^exp(XB - meanXB)
##     S = survival rate at 10 yrs
## women and men, white
## gender M:1, F:2
## bpl = BP lowering treatment, untreated or treated: 0 or 1
## Reference:
## 2013 ACC/AHA Guideline on the Assessment of Cardiovascular Risk
## A Report of the American College of Cardiology/American Heart
## Association Task Force on Practice Guidelines
## Circulation. 2014;129[suppl 2]:S49-S73

pce <- function(gender, age, tc, hdl.c, sbp, bpl, cur.smoker, dm) {
    b.surv <- c(0.9144, 0.9665)  ## baseline survival in US cohort, for white men & women
    names(b.surv) <- c("male", "female")

    coef.matrix <- array(
        dim = c(2, 11, 2),
        dimnames = list(
            c("bp.untreated", "bp.treated"),
            c(),
            c("male", "female")
            )
        )
    coef.matrix["bp.untreated", , "male"] <-
        c(12.344, 0, 11.853, -2.664, -7.990, 1.769, 1.764, 0, 7.837, -1.795, 0.658)
    coef.matrix["bp.treated", , "male"] <-
        c(12.344, 0, 11.853, -2.664, -7.990, 1.769, 1.797, 0, 7.837, -1.795, 0.658)
    coef.matrix["bp.untreated", , "female"] <-
        c(-29.799, 4.884, 13.54, -3.114, -13.578, 3.149, 1.957, 0, 7.574, -1.665, 0.661)
    coef.matrix["bp.treated", , "female"] <-
        c(-29.799, 4.884, 13.54, -3.114, -13.578, 3.149, 2.019, 0, 7.574, -1.665, 0.661)

    d.xb <- function(.df, c) {
        .df$cvs <- rowSums(data.frame(t(apply(.df, 1, function(x) {x * c}))))
        return(.df$cvs - mean(.df$cvs, na.rm = TRUE))
    }

    ar10 <- function(s, x) {
        return(1 - s**exp(x))
    }

    df <- data.frame(log(age), log(age)**2,
                     log(tc), log(age) * log(tc),
                     log(hdl.c), log(age) * log(hdl.c),
                     log(sbp), log(age)*log(sbp),
                     cur.smoker,
                     log(age) * cur.smoker,
                     dm,
                     gender, bpl)

    for (i in 1:2) {
        for (j in 1:2) {
            df[gender == i & bpl == j - 1, "risk.ascvd.10"] <-
                ar10(b.surv[i], d.xb(df[gender == i & bpl == j - 1, 1:11], coef.matrix[j, , i]))
        }
    }

    return(df$risk.ascvd.10)
}