PH251D Fall 2019 - Project 2

1. Reading in CSV file

cg <-read.csv('https://raw.githubusercontent.com/taragonmd/data/master/cage.csv', header = TRUE, sep = ',')

names(cg)

[1] "num_yes" "AUD"    

#Read in the data set and display a two-way continency table using the xtabs function. The column headings should be the AUD status.

xtabs(~num_yes + AUD, data = cg)

       AUD
num_yes   N   Y
      0 177  14
      1  22  12
      2  20  20
      3   5  43
      4   0  53

#install.packages("car")

library(car)

Loading required package: carData

cg$n.p <- cg$num_yes; cg$n.p

  [1] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
 [36] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
 [71] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 2
[106] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[141] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[176] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[211] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[246] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[281] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[316] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[351] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

cg$n.p[cg$num_yes >=3] <-"pos"
cg$n.p[cg$num_yes <=2] <-"neg"
table(cg$num_yes)

  0   1   2   3   4 
191  34  40  48  53 

table(cg$n.p)

neg pos 
265 101 

str(cg$n.p)

 chr [1:366] "pos" "pos" "pos" "pos" "pos" "pos" "pos" "pos" "pos" ...

cg$n.p<-as.factor(cg$n.p)

2. Recode integer variable into a factor with two levels

Recode the num_yes variable into a new factor variable called CAGE:

#cage <- c(AUD, ) Use the xtabs function to display 2 by 2 table of CAGE by AUD.

tab2 <-xtabs(~ n.p + AUD, data = cg); tab2

     AUD
n.p     N   Y
  neg 219  46
  pos   5  96

3. Calculate the sensitivity and specificity of the CAGE questionnaire

The sensitivity is P(CAGE = pos | AUD = Y), and the specificity is P(CAGE = neg | AUD = N).

spec<-224

4. Posterior probabilities as a function of prior probabilities

## plot of PPV vs prior probability (prior)
## use `sens` and `spec` from previous results
prior<-seq(from = 0, to = 1, by 0.05)
ppv <- (sens * prior)/(sens * prior + (1 - spec) *(1 - prior))

5. Bayesian network for probabilistic reasoning

library(bnlearn)

Attaching package: 'bnlearn'

The following object is masked from 'package:stats':

    sigma

curl <- 'https://raw.githubusercontent.com/taragonmd/data/master/cage2.csv'
aud <- read.csv(curl, header = TRUE)
summary(aud)

 AUD      CAGE    
 N:224   neg:265  
 Y:142   pos:101  

dag <- empty.graph(nodes=c("AUD", "CAGE")) # Build the Bayesian network
dag <- set.arc(dag, from = "AUD", to = "CAGE")
graphviz.plot(dag, layout = "circo") # Plot the BN

Loading required namespace: Rgraphviz

bn <- bn.fit(dag, data = aud)        # Fit the BN to the data

bn$AUD

  Parameters of node AUD (multinomial distribution)

Conditional probability table:
         N         Y 
0.6120219 0.3879781 

bn$CAGE

  Parameters of node CAGE (multinomial distribution)

Conditional probability table:
 
     AUD
CAGE           N          Y
  neg 0.97767857 0.32394366
  pos 0.02232143 0.67605634

# cpquery

library(Rgraphviz)

Loading required package: graph

Loading required package: BiocGenerics

Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:bnlearn':

    path, score

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    anyDuplicated, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter,
    Find, get, grep, grepl, intersect, is.unsorted, lapply, Map,
    mapply, match, mget, order, paste, pmax, pmax.int, pmin,
    pmin.int, Position, rank, rbind, Reduce, rownames, sapply,
    setdiff, sort, table, tapply, union, unique, unsplit, which,
    which.max, which.min

Attaching package: 'graph'

The following objects are masked from 'package:bnlearn':

    degree, nodes, nodes<-

Loading required package: grid

Attaching package: 'Rgraphviz'

The following object is masked from 'package:car':

    sp

# P(AUD ="Y" | CAGE = "pos")
(ppv1 <- cpquery(bn, event = (AUD == "Y"), evidence = (CAGE == "pos"),
n = 10^6))

[1] 0.9509407

bn$AUD$prob # view prior probabilities that we need to change

        N         Y 
0.6120219 0.3879781 

#> H
#> Pos Neg
#> 0.00001 0.99999

#Note - the prior probability was 0.32395. What if the prior probability of HIV infection was 0.05. How does the PPV change?

bn_newprior <- bn # make copy of BN model object
bn_newprior$AUD <- array(c(.05,.95), dim = dim(bn$AUD$prob),
dimnames = dimnames(bn$AUD$prob)) # assign new probs
bn_newprior$AUD$prob # confirm change to new probabilities

   N    Y 
0.05 0.95 

(ppv2 <- cpquery(bn_newprior, event = (AUD == "Y"), evidence = (CAGE == "pos"),n = 10^6))

[1] 0.9983133

6. Deconfounding: IPW with stratified contingency table

## put R code here

kidney <- read.csv("https://raw.githubusercontent.com/taragonmd/data/master/kidney.csv")
str(kidney)

'data.frame':   700 obs. of  3 variables:
 $ success   : Factor w/ 2 levels "N","Y": 2 2 2 2 2 2 2 2 2 2 ...
 $ treatment : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ stone_size: Factor w/ 2 levels "large","small": 2 2 2 2 2 2 2 2 2 2 ...

head(kidney)

table(kidney$success)

  N   Y 
138 562 

#IPW P
names(kidney) <- c('S', 'T', 'SS')

x.test <- xtabs(~ S + T + SS, data = kidney); x.test

, , SS = large

   T
S     A   B
  N  71  25
  Y 192  55

, , SS = small

   T
S     A   B
  N   6  36
  Y  81 234

tab_tss <- apply(x.test, c(2, 3), sum); tab_tss

   SS
T   large small
  A   263    87
  B    80   270

(jprob_Ytss <- prop.table(x.test) ['Y',,])

   SS
T        large      small
  A 0.27428571 0.11571429
  B 0.07857143 0.33428571

(cprob_tss <- prop.table(tab_tss, 2))

   SS
T       large     small
  A 0.7667638 0.2436975
  B 0.2332362 0.7563025

(jprob_Ytss_ipw <- jprob_Ytss/cprob_tss) #IPW matrix

   SS
T       large     small
  A 0.3577186 0.4748276
  B 0.3368750 0.4420000

(cprob_Ytss_ipw <- apply(jprob_Ytss_ipw, 1, sum)) # summation over ss

        A         B 
0.8325462 0.7788750 

r0 <- unname(cprob_Ytss_ipw['A']); r0 # P(R = yes | D = no)

[1] 0.8325462

r1 <- unname(cprob_Ytss_ipw['B']) ; r1# P(R = yes | D = yes)

[1] 0.778875

#### Average causal effects (ACEs)
#Average Causal Effects
# Causal risk difference
# Causal risk ratio
# Causal odds ratio
c(causal_RD=r1-r0, causal_RR=r1/r0, causal_OR=(r1/(1-r1))/(r0/(1-r0)))

  causal_RD   causal_RR   causal_OR 
-0.05367122  0.93553365  0.70846195 

7. Deconfounding: IPW with marginal structural model (MSM)

library(tableone)
library(sandwich) #for robust variance estimation
library(survey)

Loading required package: Matrix

Loading required package: survival

Attaching package: 'survival'

The following object is masked _by_ '.GlobalEnv':

    kidney

Attaching package: 'survey'

The following object is masked from 'package:graphics':

    dotchart

kid2 <-kidney # copy original data set with factors

kid2

kid2$S <- as.integer(kid2$S=="Y")
kid2$T <- as.integer(kid2$T=="B")
kid2$SS <- as.integer(kid2$SS=="small")

str(kidney)

'data.frame':   700 obs. of  3 variables:
 $ S : Factor w/ 2 levels "N","Y": 2 2 2 2 2 2 2 2 2 2 ...
 $ T : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ SS: Factor w/ 2 levels "large","small": 2 2 2 2 2 2 2 2 2 2 ...

head(kid2)

kid2

table(kid2$S)

  0   1 
138 562 

#### look at a table 1
table1 <- CreateTableOne(vars = "SS", strata = "T", 
                         data = kid2, test = FALSE)
#### include standardized mean difference (SMD)
print(table1,smd=TRUE)

                Stratified by T
                 0           1           SMD   
  n               350         350              
  SS (mean (SD)) 0.25 (0.43) 0.77 (0.42)  1.225

#### propensity score model
psmodel <- glm(T ~ SS, family = binomial(link ="logit"), data = kid2); psmodel

Call:  glm(formula = T ~ SS, family = binomial(link = "logit"), data = kid2)

Coefficients:
(Intercept)           SS  
     -1.190        2.323  

Degrees of Freedom: 699 Total (i.e. Null);  698 Residual
Null Deviance:      970.4 
Residual Deviance: 769.1    AIC: 773.1

ps <- predict(psmodel, type = "response"); ps

        1         2         3         4         5         6         7 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
        8         9        10        11        12        13        14 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
       15        16        17        18        19        20        21 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
       22        23        24        25        26        27        28 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
       29        30        31        32        33        34        35 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
       36        37        38        39        40        41        42 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
       43        44        45        46        47        48        49 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
       50        51        52        53        54        55        56 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
       57        58        59        60        61        62        63 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
       64        65        66        67        68        69        70 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
       71        72        73        74        75        76        77 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
       78        79        80        81        82        83        84 
0.7563025 0.7563025 0.7563025 0.7563025 0.2332362 0.2332362 0.2332362 
       85        86        87        88        89        90        91 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
       92        93        94        95        96        97        98 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
       99       100       101       102       103       104       105 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      106       107       108       109       110       111       112 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      113       114       115       116       117       118       119 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      120       121       122       123       124       125       126 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      127       128       129       130       131       132       133 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      134       135       136       137       138       139       140 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      141       142       143       144       145       146       147 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      148       149       150       151       152       153       154 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      155       156       157       158       159       160       161 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      162       163       164       165       166       167       168 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      169       170       171       172       173       174       175 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      176       177       178       179       180       181       182 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      183       184       185       186       187       188       189 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      190       191       192       193       194       195       196 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      197       198       199       200       201       202       203 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      204       205       206       207       208       209       210 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      211       212       213       214       215       216       217 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      218       219       220       221       222       223       224 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      225       226       227       228       229       230       231 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      232       233       234       235       236       237       238 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      239       240       241       242       243       244       245 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      246       247       248       249       250       251       252 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      253       254       255       256       257       258       259 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      260       261       262       263       264       265       266 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      267       268       269       270       271       272       273 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      274       275       276       277       278       279       280 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      281       282       283       284       285       286       287 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      288       289       290       291       292       293       294 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      295       296       297       298       299       300       301 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      302       303       304       305       306       307       308 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      309       310       311       312       313       314       315 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      316       317       318       319       320       321       322 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      323       324       325       326       327       328       329 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      330       331       332       333       334       335       336 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      337       338       339       340       341       342       343 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      344       345       346       347       348       349       350 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      351       352       353       354       355       356       357 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      358       359       360       361       362       363       364 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      365       366       367       368       369       370       371 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      372       373       374       375       376       377       378 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      379       380       381       382       383       384       385 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      386       387       388       389       390       391       392 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      393       394       395       396       397       398       399 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      400       401       402       403       404       405       406 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      407       408       409       410       411       412       413 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      414       415       416       417       418       419       420 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      421       422       423       424       425       426       427 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      428       429       430       431       432       433       434 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      435       436       437       438       439       440       441 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      442       443       444       445       446       447       448 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      449       450       451       452       453       454       455 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      456       457       458       459       460       461       462 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      463       464       465       466       467       468       469 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      470       471       472       473       474       475       476 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      477       478       479       480       481       482       483 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      484       485       486       487       488       489       490 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      491       492       493       494       495       496       497 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      498       499       500       501       502       503       504 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      505       506       507       508       509       510       511 
0.7563025 0.7563025 0.7563025 0.2332362 0.2332362 0.2332362 0.2332362 
      512       513       514       515       516       517       518 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      519       520       521       522       523       524       525 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      526       527       528       529       530       531       532 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      533       534       535       536       537       538       539 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      540       541       542       543       544       545       546 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      547       548       549       550       551       552       553 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      554       555       556       557       558       559       560 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      561       562       563       564       565       566       567 
0.2332362 0.2332362 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      568       569       570       571       572       573       574 
0.7563025 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      575       576       577       578       579       580       581 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      582       583       584       585       586       587       588 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      589       590       591       592       593       594       595 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      596       597       598       599       600       601       602 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      603       604       605       606       607       608       609 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      610       611       612       613       614       615       616 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      617       618       619       620       621       622       623 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      624       625       626       627       628       629       630 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      631       632       633       634       635       636       637 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      638       639       640       641       642       643       644 
0.2332362 0.2332362 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      645       646       647       648       649       650       651 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      652       653       654       655       656       657       658 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      659       660       661       662       663       664       665 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      666       667       668       669       670       671       672 
0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 0.7563025 
      673       674       675       676       677       678       679 
0.7563025 0.7563025 0.7563025 0.2332362 0.2332362 0.2332362 0.2332362 
      680       681       682       683       684       685       686 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      687       688       689       690       691       692       693 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 
      694       695       696       697       698       699       700 
0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 0.2332362 

#### create IP weights
kid2$wt <- ifelse(kid2$T==1, 1/(ps), 1/(1-ps)); kid2$wt

  [1] 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448
  [8] 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448
 [15] 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448
 [22] 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448
 [29] 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448
 [36] 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448
 [43] 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448
 [50] 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448
 [57] 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448
 [64] 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448
 [71] 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448 4.103448
 [78] 4.103448 4.103448 4.103448 4.103448 1.304183 1.304183 1.304183
 [85] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
 [92] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
 [99] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[106] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[113] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[120] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[127] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[134] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[141] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[148] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[155] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[162] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[169] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[176] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[183] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[190] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[197] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[204] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[211] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[218] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[225] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[232] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[239] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[246] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[253] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[260] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[267] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[274] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[281] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[288] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[295] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[302] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[309] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[316] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[323] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[330] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[337] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[344] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[351] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[358] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[365] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[372] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[379] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[386] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[393] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[400] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[407] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[414] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[421] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[428] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[435] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[442] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[449] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[456] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[463] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[470] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[477] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[484] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[491] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[498] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[505] 1.322222 1.322222 1.322222 4.287500 4.287500 4.287500 4.287500
[512] 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500
[519] 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500
[526] 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500
[533] 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500
[540] 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500
[547] 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500
[554] 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500
[561] 4.287500 4.287500 4.103448 4.103448 4.103448 4.103448 4.103448
[568] 4.103448 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[575] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[582] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[589] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[596] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[603] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[610] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[617] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[624] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[631] 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183 1.304183
[638] 1.304183 1.304183 1.322222 1.322222 1.322222 1.322222 1.322222
[645] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[652] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[659] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[666] 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222 1.322222
[673] 1.322222 1.322222 1.322222 4.287500 4.287500 4.287500 4.287500
[680] 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500
[687] 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500
[694] 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500 4.287500

#### apply weights to data
weighteddata <- svydesign(ids = ~ 1, data =kid2, weights = ~ wt); weighteddata

Independent Sampling design (with replacement)
svydesign(ids = ~1, data = kid2, weights = ~wt)

#### weighted table 1
weightedtable <- svyCreateTableOne(vars = "SS", strata = "T",
data = weighteddata, test = FALSE)
#### Show table with SMD
print(weightedtable, smd = TRUE)

                Stratified by T
                 0             1             SMD   
  n              700.00        700.00              
  SS (mean (SD))   0.51 (0.50)   0.51 (0.50) <0.001

#### Causal Risk Difference with weigthed GLM
glm_obj <- glm(S ~ T, weights = wt, family =
quasibinomial(link="identity"), data = kid2)
beta_ipw <- unname(coef(glm_obj))
SE <- unname(sqrt(diag(vcovHC(glm_obj, type="HC0"))))

#### get point estimate and CI for risk difference
cRD <- beta_ipw[2]
conf_level = 0.95
Z <- qnorm((1 + conf_level)/2)
cint <- beta_ipw[2] + c(-1,1) * Z * SE[2]; cint

[1] -0.12154033  0.01419789

list(causal_RD = cRD, conf_int = cint, conf_level = conf_level )

$causal_RD
[1] -0.05367122

$conf_int
[1] -0.12154033  0.01419789

$conf_level
[1] 0.95

#### Causal Risk Ratio with weighted GLM
glm_obj <- glm(S ~ T, weights = wt, family =
quasibinomial(link = log), data = kid2)
beta_ipw <- unname(coef(glm_obj))
#### To account for weighting, use asymptotic (sandwich) variance
SE <- unname(sqrt(diag(vcovHC(glm_obj, type="HC0"))))

#### get point estimate and CI for risk ratio
cRR <- exp(beta_ipw[2])
conf_level = 0.95
Z <- qnorm((1 + conf_level)/2)
cint <- exp(beta_ipw[2] + c(-1,1) * Z * SE[2])
list(causal_RR = cRR, conf_int = cint, conf_level = conf_level)

$causal_RR
[1] 0.9355336

$conf_int
[1] 0.8590792 1.0187923

$conf_level
[1] 0.95

#### Causal Odds Ratio with weighted GLM
glm_obj <- glm(S ~ T, weights = wt, family =
quasibinomial(link = logit), data = kid2)
beta_ipw <- unname(coef(glm_obj))
#### To account for weighting, use asymptotic (sandwich) variance
SE <- unname(sqrt(diag(vcovHC(glm_obj, type="HC0"))))
#### get point estimate and CI for odds ratio
cOR <- exp(beta_ipw[2])
conf_level = 0.95
Z <- qnorm((1 + conf_level)/2)
cint <- exp(beta_ipw[2] + c(-1,1) * Z * SE[2])
list(causal_OR = cOR, conf_int = cint, conf_level = conf_level)

$causal_OR
[1] 0.7084619

$conf_int
[1] 0.4617441 1.0870053

$conf_level
[1] 0.95