##############################################
#
# HOW TO USE THIS FILE
# A text version of this file (ex5.R) is also available which can just be sourced
# into the R program
# This file is intended for you to use interactively by pasting into R
# Comments are in black - red bits are commands that can be pasted into R
#
#
######### ANY TEXT AFTER THE # CHARACTER ARE TREATED AS COMMENTS
#
# data stored in c:/GR/aprojects/peas/web/exemp3/data/ex3.Rdata
# this file is c:/GR/aprojects/peas/web/exemp3/program_code/ex3_R.html
#
##############################################
#
# data are in a data frame ex5
# to get the names of the variables
#
names(ex5)
#
# the general health data is in a factor called ex5$GENHELF
# to get a table of this variable do
#
table(ex5$GENHELF)
#
# now it is time to try using the survey package
# to give R access to the survey functions do
#
library(survey)
#
# The first step is to create a SURVEY DESIGN OBJECT that
# will contain both the data and the information about the design
#
# This design has WEIGHTS (WEIGHT), but no clustering or startification
# since it is not clustered the PSU is just the individual respondent
# To represent this to the design, we set the id variable to 1
#
ex5.des<- svydesign(id=~1, weights=~WEIGHT,data=ex5)
#
# see the html help file for the function svydesign for more explanation
# the ~ sign refers to it being one of the columns of the data
# the R object now contains all the information about the design
# including the data
# To get details of the design
summary(ex5.des)
#
# now we can use the svymean commands to get first the weighted proportions
# for categorical variables
svymean(~GENHELF,ex5.des,deff=T)
svymean(~Q85A,ex5.des,deff=T)
#
# this fails for missing values so need to get subset of design without missing values
#
svymean(~Q85A,subset(ex5.des,!is.na(Q85A)),deff=T)
#
# and expressing this as percentages
#
round(print(svymean(~Q85A,subset(ex5.des,!is.na(Q85A))))*100,1)
# now the means of some continuous variables
#
svymean(~as.numeric(GENHELF),ex5.des,deff=T) # this treats it as a score not a category
svymean(~SINC,ex5.des,deff=T)
svymean(~SACC,ex5.des,deff=T)
#
# Now we recode the cannabis and amphetamine scores
# and add them into the design
CANSCORE<-as.numeric(ex5$Q85A)
CANSCORE[CANSCORE<=2]<-0
CANSCORE[CANSCORE %in% c(3,4,5)]<-1
CANSCORE[CANSCORE==6]<-0.5
ex5.des<-update(ex5.des,CANSCORE=CANSCORE)
rm(CANSCORE) # this removes vector in main workspace
#
#
AMPSCORE<-as.numeric(ex5$Q85A)
AMPSCORE[AMPSCORE<=2]<-0
AMPSCORE[AMPSCORE %in% c(3,4,5)]<-1
AMPSCORE[AMPSCORE==6]<-0.5
ex5.des<-update(ex5.des,AMPSCORE=AMPSCORE)
rm(AMPSCORE) # this removes vector in main workspace
#
# now run some regressions
#
summary(svyglm(as.numeric(GENHELF)~CANSCORE,ex5.des))
summary(svyglm(as.numeric(GENHELF)~AMPSCORE,ex5.des))
summary(svyglm(as.numeric(GENHELF)~CANSCORE+AMPSCORE,ex5.des))
#
# Now make a table to understand results
#
can.table<-svytable(~GENHELF+Q85A,ex5.des)
can.tprint<-round(sweep(can.table,2,apply(can.table,2,sum),"/")*100) # calculating column percents
can.tprint<-rbind(can.tprint,table(ex5$Q85A)) # add a row below for base numbers
#
# at any stage type tyhe name of an object to see what it is
#
dimnames(can.tprint)[[1]][6]<-'BASE'
can.tprint