Reading data into R.

Data sets can be loaded into R various ways. The simplist, is via a spreadsheet, such as a comma separated volume (csv). Using the ‘read.csv’ function. However, a working directory needs to be set-up.

# set working directly
setwd('~/OneDrive - University of Wollongong/Data/')
## setwd("C:/Users/25105517/OneDrive - University of Wollongong/Data/")

# read dataset into R, naming the data set dta.
framingham <- read.csv(file = 'framingham.csv', header = T, sep = ',')
# check the names of variable.
names(framingham)

##  [1] "randid"      "sex"         "totchol"     "age"         "sysbp"      
##  [6] "diabp"       "cursmoke"    "cigpday"     "bmi"         "diabetes"   
## [11] "bpmeds"      "heartrte"    "glucose"     "educ"        "prevchd"    
## [16] "prevap"      "prevmi"      "prevstrk"    "prevhyp"     "time"       
## [21] "period"      "hdlc"        "ldlc"        "death"       "angina"     
## [26] "hospmi"      "mi_fchd"     "anychd"      "stroke"      "cvd"        
## [31] "hyperten"    "timeap"      "timemi"      "timemifc"    "timechd"    
## [36] "timestrk"    "timecvd"     "timedth"     "timehyp"     "start.date" 
## [41] "stroke.date" "death.date"

# Look at the first ten rows.
# head(framingham)

# Look at the last ten rows.
# tail(framingham)

# View help file for dataset
# install riskCommunicator package which has the Framingham dataset.
## install.packages('riskCommunicator')
library(riskCommunicator)
## help(framingham)

# Create age groups with labels
framingham$age.group <- cut(framingham$age, breaks = c(30,39,49,59,70),
                                            labels = c('30-39','40-49',
                                                        '50-59','60-69'))
# Check lower and upper values of groups
tapply(framingham$age, framingham$age.group, min)

## 30-39 40-49 50-59 60-69 
##    32    40    50    60

tapply(framingham$age, framingham$age.group, max)

## 30-39 40-49 50-59 60-69 
##    39    49    59    70

# Convert md/dL to mmols (1 mg/dL = 18 mmols/L)
framingham$glucose.mmols <- framingham$glucose/18
framingham$totchol.mmols <- framingham$totchol/18

# Label factors with numeric values
table(framingham$sex)

## 
##    1    2 
## 5022 6605

framingham$sex <- factor(framingham$sex, level = c(1,2), 
                                         labels = c('Male','Female'))
table(framingham$sex)

## 
##   Male Female 
##   5022   6605

table(framingham$diabetes)

## 
##     0     1 
## 11097   530

framingham$diabetes <- factor(framingham$diabetes, level = c(0,1), 
                         labels = c('No','Yes'))
table(framingham$diabetes)

## 
##    No   Yes 
## 11097   530

Our specific hypothesis is related to confirming the relationship between blood pressure and the risk of stroke. Therefore, we will subset the data set.

# remove participants with a prior stroke or tx 
# for BP and only use first visit data.
names(framingham)

##  [1] "randid"        "sex"           "totchol"       "age"          
##  [5] "sysbp"         "diabp"         "cursmoke"      "cigpday"      
##  [9] "bmi"           "diabetes"      "bpmeds"        "heartrte"     
## [13] "glucose"       "educ"          "prevchd"       "prevap"       
## [17] "prevmi"        "prevstrk"      "prevhyp"       "time"         
## [21] "period"        "hdlc"          "ldlc"          "death"        
## [25] "angina"        "hospmi"        "mi_fchd"       "anychd"       
## [29] "stroke"        "cvd"           "hyperten"      "timeap"       
## [33] "timemi"        "timemifc"      "timechd"       "timestrk"     
## [37] "timecvd"       "timedth"       "timehyp"       "start.date"   
## [41] "stroke.date"   "death.date"    "age.group"     "glucose.mmols"
## [45] "totchol.mmols"

d <- subset(framingham,  bpmeds == 0 &
                         prevstrk == 0 &
                         period == 1, select = c(randid,age,totchol.mmols,
                                                 diabp,glucose.mmols,bmi,
                                                 sex,sysbp,cursmoke,educ,
                                                 cigpday,diabetes,stroke,
                                                 timestrk,death,timedth,
                                                 prevhyp,age.group))
# add n = 1, for counting
d$n <- 1
# create follow-up in years
d$timestrk.yrs <- d$timestrk/365.25
d$stroke.yes <- factor(d$stroke, levels= 0:1, labels = c('No','Yes'))
table(d$stroke.yes,d$stroke)

##      
##          0    1
##   No  3863    0
##   Yes    0  343

d$HT.yes <- factor(d$prevhyp, levels= 0:1, labels = c('No','Yes'))
table(d$HT.yes,d$prevhyp)

##      
##          0    1
##   No  2963    0
##   Yes    0 1243

# variables labels and unit
library(Hmisc)
d <- upData(d, labels = c(age = 'Age',
                          totchol.mmols = 'Total cholesterol',
                          diabp = 'Diasolic BP',
                          glucose.mmols = 'BGL',
                          sex = 'Sex',
                          sysbp = 'Systolic BP',
                          cursmoke = 'Current smoker',
                          educ = 'Education level',
                          stroke = 'Stroke',
                          timestrk = 'Follow-up',
                          timestrk.yrs = 'Follow-up'),
             units = c(age = 'yrs',
                       totchol.mmols = 'mmol/L',
                       glucose.mmols = 'mmol/L',
                       diabp = 'mmHg',
                       sysbp = 'mmHg',        
                       timestrk = 'days',
                       timestrk.yrs = 'yrs'))

## Input object size:    510712 bytes;   22 variables    4206 observations
## New object size: 501024 bytes;   22 variables    4206 observations

d$educ <- factor(d$educ, levels = 1:4, labels = c('4th grade or less',
                                                  '5th, 6th, or 7th grades',
                                                  'grade school graduate no high school',
                                                  'high school, did not graduate'))


## The Hmisc package will create Table 1. Including continuous and discrete variables. Frequencies are presented as "% (n)" and ## continuous variables as "p25/p50/p75 mean +/- SD".  

# Summary table using Hmisc package
# saving summary table a csv file

table.1 <- print(summary(sex ~ age +
                               sysbp +
                               bmi +
                               glucose.mmols +
                               diabetes +
                               HT.yes +
                               cursmoke +
                               educ +
                               timestrk.yrs,
                               method = 'reverse',
                               overall = F,
                               test = T,
                               data = d), prmsd = T, 
                                          digits = 3,
                                          prtest = 'P')

## 
## 
## Descriptive Statistics by Sex
## 
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |                                        |N   |Male                            |Female                          |P-value            |
## |                                        |    |(N=1868)                        |(N=2338)                        |                   |
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |Age [yrs]                               |4206|42.00/49.00/57.00  49.61+/- 8.70|42.00/49.00/56.00  49.62+/- 8.53|              0.808|
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |Systolic BP [mmHg]                      |4206|118.0/128.0/141.0  130.9+/- 18.8|115.5/127.5/144.0  132.3+/- 23.2|              0.894|
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |bmi                                     |4190|23.97/26.05/28.30  26.16+/- 3.40|22.53/24.71/27.64  25.47+/- 4.47|             <0.001|
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |BGL [mmol/L]                            |3822|   3.94/4.33/4.83  4.58+/-1.39  |   4.00/4.33/4.78  4.53+/-1.23  |              0.817|
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |diabetes : Yes                          |4206|            3%  (  56)          |            2%  (  51)          |              0.095|
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |HT.yes : Yes                            |4206|           31%  ( 570)          |           29%  ( 673)          |              0.222|
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |Current smoker                          |4206|           61%  (1135)          |           41%  ( 957)          |             <0.001|
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |educ : 4th grade or less                |4096|           44%  (805)           |           40%  (909)           |             <0.001|
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |    5th, 6th, or 7th grades             |    |           27%  (497)           |           32%  (724)           |                   |
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |    grade school graduate no high school|    |           13%  (228)           |           20%  (453)           |                   |
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |    high school, did not graduate       |    |           16%  (282)           |            9%  (198)           |                   |
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+
## |Follow-up [yrs]                         |4206|15.99/24.00/24.00  19.49+/- 6.59|21.09/24.00/24.00  21.02+/- 5.77|             <0.001|
## +----------------------------------------+----+--------------------------------+--------------------------------+-------------------+

write.csv(table.1, file = 'TableOne.csv', row.names = F)

Null Hypothesis Significance Testing.

\(H_0\) Mean systolic BP among non-stroke \(=\) Mean systolic BP among stroke.

\(H_i\) Mean systolic BP among non-stroke \(\neq\) Mean systolic BP among stroke.

Then a critical value is set, say from a z-statistic (i.e 1.96, < 2.5% tail)

# SHADE TAILS OF n-DISTRIBUTED PLOTS
z <- pretty(c(-3,3), 100)
ht <- dnorm(z)
data <- data.frame(z = z, ht = ht)
zc <- 1.96
t <- subset(data, z > zc)
plot(data, type = "l", 
           axes = T,
           xlab = '',
           ylab = '')

lines(data)
polygon(c(rev(t$z), t$z),c(rep(0, nrow(t)), t$ht), col = 'grey', 
                                                      border = NA)
polygon(c(rev(t$z*-1), t$z*-1), c(rep(0, nrow(t)), t$ht), 
                                        col = 'grey', border = NA)
                                        

abline(v = 0, lty = 2)
abline(v = c(-1,1), lty = 2, col = 'grey') 
text(-0.75, 0.02, '-1 SD')
text(0.75, 0.02, '+1 SD')

text(-1.75, 0.02, '-1.96 SD')
text(1.75, 0.02, '+1.96 SD')

mean(d$sysbp[d$stroke == 1]) - mean(d$sysbp[d$stroke == 0])

## [1] 13.11468

t.test(sysbp ~ stroke, data = d)

## 
##  Welch Two Sample t-test
## 
## data:  sysbp by stroke
## t = -9.1205, df = 381.19, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -15.94196 -10.28741
## sample estimates:
## mean in group 0 mean in group 1 
##        130.5967        143.7114

# EXAMPLE OF T-TEST
set.seed(12345)
n <- 3000 # large number
a <- rnorm(n, mean = 61.8, sd = 10)
b <- rnorm(n, mean = 62.7, sd = 9.5)


plot(density(a), lty = 1, col = 'blue',
     frame = F,
     main = '',
     xlab = '',
     ylab = '',
     ylim = c(0,0.075),
     lwd = 2,
     las = 1)
lines(density(b), lty = 1, col = 'black', lwd = 2)

legend(20,0.065, c('PCI - mean age 61.8 yrs (SD 10)',
                   'CABG - mean age 62.7 yrs (SD 9.5)'), 
       cex = 0.75,
       bty = 'n',
       lty = 1,
       lwd = 2,
       col = c('blue','black'))

TTest <- t.test(a,b)
TTest

## 
##  Welch Two Sample t-test
## 
## data:  a and b
## t = -3.4118, df = 5986.2, p-value = 0.0006496
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.3417023 -0.3625029
## sample estimates:
## mean of x mean of y 
##  61.76762  62.61972

str(TTest)

## List of 10
##  $ statistic  : Named num -3.41
##   ..- attr(*, "names")= chr "t"
##  $ parameter  : Named num 5986
##   ..- attr(*, "names")= chr "df"
##  $ p.value    : num 0.00065
##  $ conf.int   : num [1:2] -1.342 -0.363
##   ..- attr(*, "conf.level")= num 0.95
##  $ estimate   : Named num [1:2] 61.8 62.6
##   ..- attr(*, "names")= chr [1:2] "mean of x" "mean of y"
##  $ null.value : Named num 0
##   ..- attr(*, "names")= chr "difference in means"
##  $ stderr     : num 0.25
##  $ alternative: chr "two.sided"
##  $ method     : chr "Welch Two Sample t-test"
##  $ data.name  : chr "a and b"
##  - attr(*, "class")= chr "htest"

text(90, 0.04, paste('p-value = ', round(TTest$p.value,3), sep = ''))

n <- 30 # small number
a <- rnorm(n, mean = 61.8, sd = 10)
b <- rnorm(n, mean = 62.7, sd = 9.5)


plot(density(a), lty = 1, col = 'blue',
     frame = F,
     main = '',
     xlab = '',
     ylab = '',
     ylim = c(0,0.075),
     lwd = 2,
     las = 1)
lines(density(b), lty = 1, col = 'black', lwd = 2)

legend(50,0.065, c('PCI - mean age 61.8 yrs (SD 10)',
                   'CABG - mean age 62.7 yrs (SD 9.5)'), 
       cex = 0.75,
       bty = 'n',
       lty = 1,
       lwd = 2,
       col = c('blue','black'))

TTest <- t.test(a,b)
TTest

## 
##  Welch Two Sample t-test
## 
## data:  a and b
## t = -0.28393, df = 54.736, p-value = 0.7775
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.455870  4.101878
## sample estimates:
## mean of x mean of y 
##  61.79375  62.47075

str(TTest)

## List of 10
##  $ statistic  : Named num -0.284
##   ..- attr(*, "names")= chr "t"
##  $ parameter  : Named num 54.7
##   ..- attr(*, "names")= chr "df"
##  $ p.value    : num 0.778
##  $ conf.int   : num [1:2] -5.46 4.1
##   ..- attr(*, "conf.level")= num 0.95
##  $ estimate   : Named num [1:2] 61.8 62.5
##   ..- attr(*, "names")= chr [1:2] "mean of x" "mean of y"
##  $ null.value : Named num 0
##   ..- attr(*, "names")= chr "difference in means"
##  $ stderr     : num 2.38
##  $ alternative: chr "two.sided"
##  $ method     : chr "Welch Two Sample t-test"
##  $ data.name  : chr "a and b"
##  - attr(*, "class")= chr "htest"

text(80, 0.04, paste('p-value = ', round(TTest$p.value,3), sep = ''))

## Type I and Type II error

NEJM 2000 ARDSNet lower tidal volume (6 mL/kg) versus traditional 12mL/kg. The primary outome death before discharge and breathing without assitance. Results were tx group 134/432 (31.0%) versus control 171/429 (39.8%).

## setwd('~/OneDrive - University of Wollongong/Data/')
## setwd("C:/Users/25105517/OneDrive - University of Wollongong/Data/")
source("~/OneDrive - University of Wollongong/Data/epiTablesV3.R")

## 
## ------------------------------------------------------------------------------
## Epi Tables
## R epiTab Function
## Input values epiTab(a,b,c,d)
## a,b,c,d values from two-by-two table (rows exposure and columns outcome)
## -------------------------------------------------------------------------------
## ------------------------------------------------------------------------
## Trial Data:
##             Outcome
## Intervention  Yes   No Total
##        Yes    482 1622  2104
##        No    1110 3211  4321
##        Total 1592 4833  6425
## 
## Rate of outcome among intervention group (p1) = 0.2291 (95% CI: 21.11, 24.7)
## 
## Rate of outcome among control group (p2) = 0.2569 (95% CI: 24.39, 26.99)
## 
## Risk Difference = -2.78 (95% CI: -5, -0.56)
## 
## Hazard Risk = 0.876 95% CI = 0.787 to 0.975
## Rate Ratio = 0.892 95% CI = 0.812 to 0.979
## Odds Ratio = 0.86 95% CI = 0.761 to 0.972
## chi2 p-value = 0.0154
## Fisher's exact p-value = 0.0163
## ------------------------------------------------------------------------
## 
## ------------------------------------------------------------------------
## Trial Data:
##             Outcome
## Intervention Yes  No Total
##        Yes    95 229   324
##        No    283 400   683
##        Total 378 629  1007
## 
## Rate of outcome among intervention group (p1) = 0.2932 (95% CI: 24.34, 34.3)
## 
## Rate of outcome among control group (p2) = 0.4143 (95% CI: 37.73, 45.14)
## 
## Risk Difference = -12.11 (95% CI: -18.3, -5.92)
## 
## Hazard Risk = 0.649 95% CI = 0.514 to 0.819
## Rate Ratio = 0.708 95% CI = 0.585 to 0.857
## Odds Ratio = 0.586 95% CI = 0.441 to 0.778
## chi2 p-value = 2e-04
## Fisher's exact p-value = 2e-04
## ------------------------------------------------------------------------

epiTab(134,298,171,258)

## 
## ------------------------------------------------------------------------
## Trial Data:
##             Outcome
## Intervention Yes  No Total
##        Yes   134 298   432
##        No    171 258   429
##        Total 305 556   861
## 
## Rate of outcome among intervention group (p1) = 0.3102 (95% CI: 26.64, 35.39)
## 
## Rate of outcome among control group (p2) = 0.3986 (95% CI: 35.21, 44.51)
## 
## Risk Difference = -8.84 (95% CI: -15.21, -2.47)
## 
## Hazard Risk = 0.73 95% CI = 0.582 to 0.915
## Rate Ratio = 0.778 95% CI = 0.648 to 0.934
## Odds Ratio = 0.678 95% CI = 0.512 to 0.898
## chi2 p-value = 0.0067
## Fisher's exact p-value = 0.0069
## ------------------------------------------------------------------------

An example of potentially Type II error - NEJM 2018 EOLIA ECMO for ARDS primary outome 60-day mortality tx 44/124 (35%) control 57/125 (46%)

epiTab(44,80,57,68)

## 
## ------------------------------------------------------------------------
## Trial Data:
##             Outcome
## Intervention Yes  No Total
##        Yes    44  80   124
##        No     57  68   125
##        Total 101 148   249
## 
## Rate of outcome among intervention group (p1) = 0.3548 (95% CI: 26.98, 43.99)
## 
## Rate of outcome among control group (p2) = 0.456 (95% CI: 36.78, 54.42)
## 
## Risk Difference = -10.12 (95% CI: -22.31, 2.07)
## 
## Hazard Risk = 0.72 95% CI = 0.486 to 1.067
## Rate Ratio = 0.778 95% CI = 0.574 to 1.055
## Odds Ratio = 0.656 95% CI = 0.394 to 1.091
## chi2 p-value = 0.1041
## Fisher's exact p-value = 0.1217
## ------------------------------------------------------------------------

# original sample size calculation was based on 40% vs 60%
# 60-day mortality
power.prop.test(p1 = 0.4, p2 = 0.6, sig.level = 0.05, power = 0.8)

## 
##      Two-sample comparison of proportions power calculation 
## 
##               n = 96.92364
##              p1 = 0.4
##              p2 = 0.6
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

p1 <- 0.6
q1 <- 1 - p1
p2 <- 0.4
q2 <- 1 - p2
p.bar <- (p1 + p2)/2
q.bar <- 1 - p.bar

# Sample size
n <- 2*((p.bar * q.bar) * (1.96 + 0.84)^2)  / (p1 - p2)^2
n

## [1] 98

power.prop.test(p1 = 0.36, p2 = 0.46, sig.level = 0.05, power = 0.8)

## 
##      Two-sample comparison of proportions power calculation 
## 
##               n = 378.5477
##              p1 = 0.36
##              p2 = 0.46
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

power.prop.test(p1 = 0.36, p2 = 0.46, sig.level = 0.05, n = 125) # power

## 
##      Two-sample comparison of proportions power calculation 
## 
##               n = 125
##              p1 = 0.36
##              p2 = 0.46
##       sig.level = 0.05
##           power = 0.3615169
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

# Age and stroke
boxplot(age ~ stroke.yes, data = d,
        col = c('blue','red'),
        xlab = 'Stroke',
        ylab = 'Age (years)',
        las = 1)

tapply(d$age, d$stroke, mean)

##        0        1 
## 49.13461 55.07289

tapply(d$age, d$stroke, sd)

##        0        1 
## 8.504941 7.768804

# Age and Systolic BP
plot(d$age, d$sysbp, pch = 16,
                     xlab = 'Age (yrs)',
                     ylab = 'Systolic BP (mmHg)',
                     main = '')

plot(d$age[d$stroke == 0], d$sysbp[d$stroke == 0], pch = 16,
                     col = 'blue',
                     xlab = 'Age (yrs)',
                     ylab = 'Systolic BP (mmHg)',
                     main = '',
                     ylim = c(80,300),
                     las = 1)
points(d$age[d$stroke == 1], d$sysbp[d$stroke == 1], pch = 16,
                             col = 'red')
abline(lsfit(d$age,d$sysbp), lty = 2, lwd = 2)

# Cigs per day and stroke
boxplot(cigpday ~ stroke.yes, data = d,
        col = c('blue','red'),
        xlab = 'Stroke',
        ylab = 'Age (years)',
        las = 1)

tapply(d$cigpday, d$stroke, mean, na.rm = T)

##        0        1 
## 9.105085 9.170088

tapply(d$cigpday, d$stroke, sd, na.rm = T)

##        0        1 
## 12.00131 12.06747

# table of rates of stroke based 
# systolic HT
library(Epi)
d$sysHT <- ifelse(d$sysbp > 140, 1, 0)
d$sysHT <- factor(d$sysHT, levels = 0:1, labels = c('No','Yes'))
print(stat.table(list(sysHT), contents = list(N = count(),
                                        'Stroke (n)' = sum(stroke),
                                        'Stroke (%)' = ratio(stroke,n,100),
                                        'Stroke (rate)' = ratio(stroke,timestrk,365.25),
                                        'follow-up (Total years)' = sum(timestrk.yrs)),
                                         data = d,
                                         margin = T),
                                         digits = 3)

##  --------------------------------------------------- 
##  sysHT          N  Stroke  Stroke  Stroke follow-up  
##                       (n)     (%)  (rate)    (Total  
##                                              years)  
##  --------------------------------------------------- 
##  No      3035.000 178.000   5.865   0.003 64518.680  
##  Yes     1171.000 165.000  14.091   0.008 21014.623  
##                                                      
##  Total   4206.000 343.000   8.155   0.004 85533.303  
##  ---------------------------------------------------

tapply(d$timestrk.yrs, d$sysHT, sum)

##       No      Yes 
## 64518.68 21014.62

# rate with 95% CI
library(epiDisplay)
stroke.tab <- data.frame(cbind(tapply(d$stroke, d$sysHT, sum), tapply(d$n, d$sysHT, sum)))
names(stroke.tab) <- c('Strokes','N')
head(stroke.tab)

##     Strokes    N
## No      178 3035
## Yes     165 1171

ci.binomial(stroke.tab$Strokes, stroke.tab$N)

##   events total probability          se exact.lower95ci exact.upper95ci
## 1    178  3035  0.05864909 0.004265079      0.05054965      0.06761068
## 2    165  1171  0.14090521 0.010167300      0.12146029      0.16216781

# incidence based on person years
stroke.tab <- data.frame(cbind(tapply(d$stroke, d$sysHT, sum), tapply(d$timestrk.yrs, d$sysHT, sum)))
names(stroke.tab) <- c('Strokes','pyrs')
head(stroke.tab)

##     Strokes     pyrs
## No      178 64518.68
## Yes     165 21014.62

ci.poisson(stroke.tab$Strokes, stroke.tab$pyrs)

##   events person.time   incidence           se exact.lower95ci exact.upper95ci
## 1    178    64518.68 0.002758891 0.0002067876     0.002367128     0.003197555
## 2    165    21014.62 0.007851676 0.0006112521     0.006697479     0.009147202

# risk of stroke
# Rate Ratio
## setwd('~/OneDrive - University of Wollongong/Data/')
## setwd("C:/Users/25105517/OneDrive - University of Wollongong/Data/")
source('~/OneDrive - University of Wollongong/Data/epiTablesV3.R')

## 
## ------------------------------------------------------------------------------
## Epi Tables
## R epiTab Function
## Input values epiTab(a,b,c,d)
## a,b,c,d values from two-by-two table (rows exposure and columns outcome)
## -------------------------------------------------------------------------------
## ------------------------------------------------------------------------
## Trial Data:
##             Outcome
## Intervention  Yes   No Total
##        Yes    482 1622  2104
##        No    1110 3211  4321
##        Total 1592 4833  6425
## 
## Rate of outcome among intervention group (p1) = 0.2291 (95% CI: 21.11, 24.7)
## 
## Rate of outcome among control group (p2) = 0.2569 (95% CI: 24.39, 26.99)
## 
## Risk Difference = -2.78 (95% CI: -5, -0.56)
## 
## Hazard Risk = 0.876 95% CI = 0.787 to 0.975
## Rate Ratio = 0.892 95% CI = 0.812 to 0.979
## Odds Ratio = 0.86 95% CI = 0.761 to 0.972
## chi2 p-value = 0.0154
## Fisher's exact p-value = 0.0163
## ------------------------------------------------------------------------
## 
## ------------------------------------------------------------------------
## Trial Data:
##             Outcome
## Intervention Yes  No Total
##        Yes    95 229   324
##        No    283 400   683
##        Total 378 629  1007
## 
## Rate of outcome among intervention group (p1) = 0.2932 (95% CI: 24.34, 34.3)
## 
## Rate of outcome among control group (p2) = 0.4143 (95% CI: 37.73, 45.14)
## 
## Risk Difference = -12.11 (95% CI: -18.3, -5.92)
## 
## Hazard Risk = 0.649 95% CI = 0.514 to 0.819
## Rate Ratio = 0.708 95% CI = 0.585 to 0.857
## Odds Ratio = 0.586 95% CI = 0.441 to 0.778
## chi2 p-value = 2e-04
## Fisher's exact p-value = 2e-04
## ------------------------------------------------------------------------

epiTab(165,1171-165,178,3035-178)

## 
## ------------------------------------------------------------------------
## Trial Data:
##             Outcome
## Intervention Yes   No Total
##        Yes   165 1006  1171
##        No    178 2857  3035
##        Total 343 3863  4206
## 
## Rate of outcome among intervention group (p1) = 0.1409 (95% CI: 12.1, 16.09)
## 
## Rate of outcome among control group (p2) = 0.0586 (95% CI: 5.03, 6.7)
## 
## Risk Difference = 8.23 (95% CI: 6.06, 10.39)
## 
## Hazard Risk = 2.513 95% CI = 2.033 to 3.106
## Rate Ratio = 2.403 95% CI = 1.966 to 2.937
## Odds Ratio = 2.633 95% CI = 2.105 to 3.293
## chi2 p-value = 0
## Fisher's exact p-value = 0
## ------------------------------------------------------------------------

(sum(d$stroke[d$sysHT == 'Yes']) / sum(d$n[d$sysHT == 'Yes'])) / 
(sum(d$stroke[d$sysHT == 'No']) / sum(d$n[d$sysHT == 'No']))

## [1] 2.402513

Confidence intervals

Confidence intervals can be found all over statistics. They provide an interval likely to include the true population parameter we’re trying to estimate, allowing us to express estimated values from sample data with some confidence.

Depending on the situation, there are numerous methods for calculating them.The following formula is used to compute it:

    $$Confidence Interval = (point estimate)+/-(critical value)*(standard error)$$

Module 4. Inferential statistics 2: Correlation and regression, adjusting for confounders.

This module will initially explore correlation and simple linear regression. We will then extend these concepts to multi-linear regression, and generalised-linear-models that are use to analysis various outcomes, such as continuous, binomial, and rates.

# Correlation and linear regression

set.seed(12345)
d.samp <-subset(d, randid %in% sample(d$randid, 150))

par(mfrow = c(1,1))
plot(sysbp ~ age, data = d.samp,
                  las = 1,
                  xlab = 'Age (yrs)',
                  ylab = 'Systolic BP (mmHg)')

summary(lm(sysbp ~ age, data = d.samp))

## 
## Call:
## lm(formula = sysbp ~ age, data = d.samp)
## 
## Residuals:
## Systolic BP [mmHg] 
##     Min      1Q  Median      3Q     Max 
## -49.120 -14.068  -1.351   9.512  91.466 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  87.9392     9.7463   9.023 8.91e-16 ***
## age           0.8622     0.1925   4.479 1.49e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.21 on 148 degrees of freedom
## Multiple R-squared:  0.1194, Adjusted R-squared:  0.1134 
## F-statistic: 20.06 on 1 and 148 DF,  p-value: 1.489e-05

abline(a = 87.9392, b = 0.8622, lty = 2)

summary(lm(sysbp ~ I(age-35), data = d.samp))

## 
## Call:
## lm(formula = sysbp ~ I(age - 35), data = d.samp)
## 
## Residuals:
## Systolic BP [mmHg] 
##     Min      1Q  Median      3Q     Max 
## -49.120 -14.068  -1.351   9.512  91.466 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 118.1164     3.3385  35.380  < 2e-16 ***
## I(age - 35)   0.8622     0.1925   4.479 1.49e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.21 on 148 degrees of freedom
## Multiple R-squared:  0.1194, Adjusted R-squared:  0.1134 
## F-statistic: 20.06 on 1 and 148 DF,  p-value: 1.489e-05

cor(d.samp$age, d.samp$sysbp)

## [1] 0.3455124

cov(d.samp$age, d.samp$sysbp) / (var(d.samp$age) * var(d.samp$sysbp))^0.5

## [1] 0.3455124

tapply(d$sysbp, d$stroke, mean)

##        0        1 
## 130.5967 143.7114

summary(lm(sysbp ~ stroke, data = d))

## 
## Call:
## lm(formula = sysbp ~ stroke, data = d)
## 
## Residuals:
## Systolic BP [mmHg] 
##     Min      1Q  Median      3Q     Max 
## -49.711 -14.597  -3.597  10.903 151.289 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 130.5967     0.3386  385.66   <2e-16 ***
## stroke       13.1147     1.1858   11.06   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.05 on 4204 degrees of freedom
## Multiple R-squared:  0.02827,    Adjusted R-squared:  0.02804 
## F-statistic: 122.3 on 1 and 4204 DF,  p-value: < 2.2e-16

summary(lm(sysbp ~ stroke - 1, data = d))

## 
## Call:
## lm(formula = sysbp ~ stroke - 1, data = d)
## 
## Residuals:
## Systolic BP [mmHg] 
##    Min     1Q Median     3Q    Max 
## -49.71 113.00 125.00 139.88 243.00 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)    
## stroke  143.711      6.854   20.97   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 126.9 on 4205 degrees of freedom
## Multiple R-squared:  0.09467,    Adjusted R-squared:  0.09445 
## F-statistic: 439.7 on 1 and 4205 DF,  p-value: < 2.2e-16

regress.display(lm(sysbp ~ stroke, data = d))

## Linear regression predicting sysbp
##  
##                 Coeff.(95%CI)        P(t-test) P(F-test)
## stroke: 1 vs 0  13.11 (10.79,15.44)  < 0.001   < 0.001  
##                                                         
## No. of observations = 4206

# crude models 
plot(stroke ~ sysbp, data = d, las = 1)

plot(sysbp ~ stroke.yes, data = d, las = 1)

stat.table(HT.yes, contents = list(D = sum(stroke),
                                   N = sum(n),
                                   R = ratio(stroke,n)),
                                   margin = T,
                                   data = d)

##  --------------------------------- 
##  HT.yes         D       N       R  
##  --------------------------------- 
##  No        157.00 2963.00    0.05  
##  Yes       186.00 1243.00    0.15  
##                                    
##  Total     343.00 4206.00    0.08  
##  ---------------------------------

# logit model
summary(glm(stroke ~ HT.yes, data = d, family = binomial(link='logit')))

## 
## Call:
## glm(formula = stroke ~ HT.yes, family = binomial(link = "logit"), 
##     data = d)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5694  -0.5694  -0.3300  -0.3300   2.4239  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.88327    0.08201  -35.16   <2e-16 ***
## HT.yesYes    1.14583    0.11423   10.03   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2376.7  on 4205  degrees of freedom
## Residual deviance: 2277.3  on 4204  degrees of freedom
## AIC: 2281.3
## 
## Number of Fisher Scoring iterations: 5

# predicted rate exp(a+bx)/(a+bx)+1
exp(-2.88327 + (1*1.1458)) / (exp(-2.88327 + (1*1.1458)) + 1)

## [1] 0.1496346

# log odds
# observed odds ratio
(0.15/0.85)/(0.05/0.95)

## [1] 3.352941

# predicted OR
exp(1.14583)

## [1] 3.145051

# Rate Ratio
tapply(d$stroke, d$HT.yes, mean)

##         No        Yes 
## 0.05298684 0.14963797

0.14963797/0.05298684

## [1] 2.824059

glm.stroke <- glm(cbind(stroke,n) ~ sysHT, data = d, family = binomial)
print(logistic.display(glm.stroke), digits = 3)

## 
## Logistic regression predicting cbind(stroke, n) 
##  
##                  OR(95%CI)     P(Wald's test) P(LR-test)
## sysHT: Yes vs No 2.4 (1.92,3)  < 0.001        < 0.001   
##                                                         
## Log-likelihood = -949.6711
## No. of observations = 4206
## AIC value = 1903.3422

glm.stroke <- glm(stroke ~ sysHT + offset(n), data = d, family = poisson)
glm.stroke <- glm(stroke ~ sysHT, offset(n), data = d, family = poisson)
print(idr.display(glm.stroke), digits = 3)

## 
## 
## Poisson regression predicting stroke 
##  
##                  IDR(95%CI)       P(Wald's test) P(LR-test)
## sysHT: Yes vs No 2.4 (1.94,2.97)  < 0.001        < 0.001   
##                                                            
## Log-likelihood = -1171.1858
## No. of observations = 4206
## AIC value = 2346.3716

cph <- coxph(Surv(timestrk, stroke) ~ I(sysbp/10), data = d)
summary(cph)

## Call:
## coxph(formula = Surv(timestrk, stroke) ~ I(sysbp/10), data = d)
## 
##   n= 4206, number of events= 343 
## 
##                coef exp(coef) se(coef)     z Pr(>|z|)    
## I(sysbp/10) 0.28278   1.32681  0.02013 14.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             exp(coef) exp(-coef) lower .95 upper .95
## I(sysbp/10)     1.327     0.7537     1.275      1.38
## 
## Concordance= 0.684  (se = 0.015 )
## Likelihood ratio test= 157  on 1 df,   p=<2e-16
## Wald test            = 197.3  on 1 df,   p=<2e-16
## Score (logrank) test = 196.6  on 1 df,   p=<2e-16

# adjusted models 
# freq data
stroke.tab <- data.frame(cbind(tapply(d$stroke, d$sysHT, sum), tapply(d$n, d$sysHT, sum)))
names(stroke.tab) <- c('Strokes','N')
ci.binomial(stroke.tab$Strokes, stroke.tab$N)

##   events total probability          se exact.lower95ci exact.upper95ci
## 1    178  3035  0.05864909 0.004265079      0.05054965      0.06761068
## 2    165  1171  0.14090521 0.010167300      0.12146029      0.16216781

stroke.tab <- data.frame(cbind(tapply(d$stroke, d$sex, sum), tapply(d$n, d$sex, sum)))
names(stroke.tab) <- c('Strokes','N')
ci.binomial(stroke.tab$Strokes, stroke.tab$N)

##   events total probability          se exact.lower95ci exact.upper95ci
## 1    167  1868  0.08940043 0.006601537      0.07683967       0.1032664
## 2    176  2338  0.07527802 0.005456541      0.06490470       0.0867278

stroke.tab <- data.frame(cbind(tapply(d$stroke, d$age.group, sum), tapply(d$n, d$age.group, sum)))
names(stroke.tab) <- c('Strokes','N')
ci.binomial(stroke.tab$Strokes, stroke.tab$N)

##   events total probability          se exact.lower95ci exact.upper95ci
## 1     11   553  0.01989150 0.005937564      0.00996962      0.03531338
## 2     68  1642  0.04141291 0.004916965      0.03229793      0.05220926
## 3    148  1307  0.11323642 0.008765147      0.09655669      0.13168263
## 4    116   704  0.16477273 0.013981658      0.13809602      0.19428352

table(d$age.group, d$stroke.yes)

##        
##          Yes   No
##   30-39   11  542
##   40-49   68 1574
##   50-59  148 1159
##   60-69  116  588

round(epi.conf(dat = table(d$age.group, d$stroke.yes), ctype = "prop.single"), digits = 3)

##         est    se lower upper
## 30-39 0.020 0.006 0.011 0.035
## 40-49 0.041 0.005 0.033 0.052
## 50-59 0.113 0.009 0.097 0.132
## 60-69 0.165 0.014 0.139 0.194

d$sex <- relevel(d$sex, ref = 'Female')
glm.stroke <- glm(cbind(stroke,n) ~ sysHT + sex + age.group, data = d, family = binomial)
print(logistic.display(glm.stroke), digits = 3)

## 
## Logistic regression predicting cbind(stroke, n) 
##  
##                       crude OR(95%CI)    adj. OR(95%CI)     P(Wald's test)
## sysHT: Yes vs No      2.4 (1.92,3)       1.72 (1.36,2.17)   < 0.001       
##                                                                           
## sex: Male vs Female   1.19 (0.95,1.48)   1.27 (1.01,1.59)   0.038         
##                                                                           
## age.group: ref.=30-39                                                     
##    40-49              2.08 (1.09,3.96)   1.96 (1.03,3.74)   0.041         
##    50-59              5.69 (3.06,10.59)  4.84 (2.59,9.05)   < 0.001       
##    60-69              8.28 (4.42,15.53)  6.49 (3.42,12.28)  < 0.001       
##                                                                           
##                       P(LR-test)
## sysHT: Yes vs No      < 0.001   
##                                 
## sex: Male vs Female   < 0.001   
##                                 
## age.group: ref.=30-39 < 0.001   
##    40-49                        
##    50-59                        
##    60-69                        
##                                 
## Log-likelihood = -903.8435
## No. of observations = 4206
## AIC value = 1819.687

# Using follow-up time
# Kaplan-Meier estimates
km.0 <- survfit(Surv(timestrk.yrs, stroke) ~ sysHT, data = d)
summary(km.0, times = c(5,10,20))

## Call: survfit(formula = Surv(timestrk.yrs, stroke) ~ sysHT, data = d)
## 
##                 sysHT=No 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5   2945      13    0.996 0.00120        0.993        0.998
##    10   2819      20    0.989 0.00194        0.985        0.993
##    20   2374      91    0.955 0.00398        0.947        0.963
## 
##                 sysHT=Yes 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5   1082      24    0.979 0.00427        0.971        0.987
##    10    950      36    0.944 0.00702        0.931        0.958
##    20    639      74    0.858 0.01154        0.836        0.881

par(mfrow = c(1,2))
plot(km.0, lty = c(1,2),
           las = 1,
           col = c('blue','red'),
           xlab = 'Follow-up (yrs)',
           ylab = 'Probability of avoiding stroke')

legend('bottomleft', c('no-HT','HT'),
                     lty = c(1,2),
                     col = c('blue','red'),
                     bty = 'n')

plot(km.0, fun = 'event',
           lty = c(1,2),
           col = c('blue','red'),
           las = 1,
           xlab = 'Follow-up (yrs)',
           ylab = 'Probability of stroke')

legend('topleft', c('no-HT','HT'),
                  lty = c(1,2),
                  col = c('blue','red'),
                  bty = 'n')

cph.0 <- coxph(Surv(timestrk.yrs, stroke) ~ sysHT, data = d)
summary(cph.0)

## Call:
## coxph(formula = Surv(timestrk.yrs, stroke) ~ sysHT, data = d)
## 
##   n= 4206, number of events= 343 
## 
##            coef exp(coef) se(coef)     z Pr(>|z|)    
## sysHTYes 1.1139    3.0464   0.1083 10.29   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## sysHTYes     3.046     0.3283     2.464     3.767
## 
## Concordance= 0.628  (se = 0.013 )
## Likelihood ratio test= 99.52  on 1 df,   p=<2e-16
## Wald test            = 105.8  on 1 df,   p=<2e-16
## Score (logrank) test = 117.1  on 1 df,   p=<2e-16

newData <- data.frame(expand.grid(stroke = 0:1,
                                  sysHT = c('Yes','No'),
                                  timestrk.yrs = c(5,10,20)))

newData

##    stroke sysHT timestrk.yrs
## 1       0   Yes            5
## 2       1   Yes            5
## 3       0    No            5
## 4       1    No            5
## 5       0   Yes           10
## 6       1   Yes           10
## 7       0    No           10
## 8       1    No           10
## 9       0   Yes           20
## 10      1   Yes           20
## 11      0    No           20
## 12      1    No           20

yhat <- predict(cph.0, newdata = newData, type = 'expected')
pred.stroke <- cbind(newData, yhat)
subset(pred.stroke, stroke == 1 & timestrk.yrs %in% c(10,20))

##    stroke sysHT timestrk.yrs       yhat
## 6       1   Yes           10 0.04631465
## 8       1    No           10 0.01520324
## 10      1   Yes           20 0.14772021
## 12      1    No           20 0.04849063

par(mfrow = c(1,1),
    mar = c(8,4,2,8),
    cex = 0.9)
plot(km.0, fun = 'event',
     lty = c(2,2),
     col = c('blue','red'),
     las = 1,
     lwd = 2,
     xlab = 'Follow-up (yrs)',
     ylab = 'Probability of stroke',
     ylim = c(0,0.15),
     las = 1,
     frame = F,
     axes = F,
     xmax = 20)
axis(1, at = seq(0,20,1), label = seq(0,20,1))
axis(2, at = seq(0,0.15,0.01), label = seq(0,0.15,0.01), las = 1)

summary(km.0, times = c(5,10,20))

## Call: survfit(formula = Surv(timestrk.yrs, stroke) ~ sysHT, data = d)
## 
##                 sysHT=No 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5   2945      13    0.996 0.00120        0.993        0.998
##    10   2819      20    0.989 0.00194        0.985        0.993
##    20   2374      91    0.955 0.00398        0.947        0.963
## 
##                 sysHT=Yes 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5   1082      24    0.979 0.00427        0.971        0.987
##    10    950      36    0.944 0.00702        0.931        0.958
##    20    639      74    0.858 0.01154        0.836        0.881

mtext(side = 4,at = 1-0.858,  col = "red", text = 'Hypertension', las = 1)
mtext(side = 4,at = 1-0.955,  col = "blue", text = 'No-hypertension', las = 1)

# Add no. at risk
n.surv <- summary(km.0, times = seq(0,20,1))
## n.surv
n.risk <- matrix(n.surv$n.risk, byrow = T, nrow = 2)
## n.risk

mtext(side = 1, at = -0.1, line = 4, text = 'no. at risk \nHT group', outer = F, cex = 0.7, font = 1, adj = 0)
mtext(side = 1, at = c(-1, seq(0,20,1)), line = 5, text = c('HT', n.risk[1,]), outer = F, cex = 0.7, adj = 0, col = 1)
mtext(side = 1, at = c(-1, seq(0,20,1)), line = 6, text = c('No HT', n.risk[2,]), outer = F, cex = 0.7, adj = 0, col = 2)

Table 1. Characteristics of study participants.

Using the Hmisc ‘summary’ function we have created the above Table, we have included test = TRUE, so the function will compare the continuous variables using (un-paired T-Test), and the frequency data using a Pearson’s chi-squared test.

Table 2. Hypertension and risk of stroke.

• adjusted \(p\)-values.

Some useful links:

BMJ Statistics at square one.

https://www.bmj.com/about-bmj/resources-readers/publications/statistics-square-one

Bland and Altman Statistical notes in the BMJ.

https://www-users.york.ac.uk/~mb55/pubs/pbstnote.htm