Edps/Psych/Soc 584
# Fall 2018               last revised:  June 6 2018
# C.J.Anderson
#
#  Illustrates computations for:
#   (a) Data input as freqency and cross-classification 
#   (b) Estimates of differences proportions, relative risk and odds rati0

library(xtable)        # To output tables in LaTex
library(Epi)           # For getting relative risk, odds ratio, difference of p
library(MASS)          # Has lots of things, below needed for loglm
library(vcd)           # Needed for assocstats
library(vcdExtra)      # For oddsratio()&  odds ratio() commands

# There are a number of different ways to enter data as tables.
# This one is simple and easy for 2 x 2 tables

# Enter in data frequency form (case form)
# -- set-up variables names and levels
var.levels <- expand.grid( cold=c("yes","no"), 
                           treatment=c("vitamin C","placebo"))

# -- create a data frame with  2 x 2 = 4 cases				   
( colds.freq <- data.frame(var.levels, count=c(17,122,31,109)) )


# Data set-up:  Put data into tabular form
( colds.tab <- xtabs(count ~ treatment + cold,data=colds.freq) )

# For LaTex output 

xtable(colds.tab)


# (b) 

# Difference of proportions, Relative-Risk and Odds Ratios
# Package Epi

# data should be in matrix form -- this keeps row & column names and level values
( colds.mat <- as.matrix(colds.tab) )

# chisquar **with continuity correction**
twoby2(colds.mat)

# chisquare without correction (this is just fine!)
loglm(count ~ cold + treatment, data=colds.freq)

# Or more 

assocstats(colds.tab)

OR <- oddsratio(colds.tab)
summary(OR)
OR                         # log(odds ratio)                         
confint(OR,level=.95)      # 95% CI for log (odds ratio)
exp(confint(OR,level=.95) ) # 95% CI for odds ratio