Vincent Schubert Malbas ®

SAS code

libname mydata "/courses/u_coursera.org1/i_1006328/c_5333" access=readonly;

DATA new; set mydata.nesarc_pds;

/*major depressive disorder ~ Sex, average daily quantity of alcohol consumed and cigarettes smoked in past 12 months, and use of sedatives, tranquilizers, cannabis, opioids, amphetamines, cocaine, heroine, hallucinogens, and inhalants*/

LABEL MAJORDEPLIFE="major depressive disorder (lifetime)"
SEX="biological sex"
S2AQ8B="alcohol consumed daily" 
S3AQ3C1="cigarettes smoked in past 12 months"
S3BQ1A1="use of sedatives"
S3BQ1A2="use of  tranquilizers"
S3BQ1A3="use of opioids"
S3BQ1A4="use of amphetamines"
S3BQ1A5="use of cannabis"
S3BQ1A6= "use of cocaine"
S3BQ1A7= "use of hallucinogens"
S3BQ1A8= "use of inhalants"
S3BQ1A9A="use of heroin";

/* DATA MUNGING */

/* for the  following variables: alcohol consumed, cigarettes smoked */
 /* change [9] to [NA] */

IF S2AQ8B=99 THEN S2AQ8B=.; /*alcohol*/
IF S3AQ3C1=99 THEN S3AQ3C1=.; /*cigarette*/

/* for the  following variables: "sex", "sedatives", "tranquilizers", "cannabis", "opioids", "amphetamines", "cocaine", "heroine", "hallucinogens", "inhalants" */
 /* change the binomial setting [1,2] --> [1,0] */
 /* change [9] to [NA] */

IF SEX=2 THEN SEX=0;

IF S3BQ1A1=2 THEN S3BQ1A1=0;
IF S3BQ1A2=2 THEN S3BQ1A2=0;
IF S3BQ1A3=2 THEN S3BQ1A3=0;
IF S3BQ1A4=2 THEN S3BQ1A4=0;
IF S3BQ1A5=2 THEN S3BQ1A5=0;
IF S3BQ1A6=2 THEN S3BQ1A6=0;
IF S3BQ1A7=2 THEN S3BQ1A7=0;
IF S3BQ1A8=2 THEN S3BQ1A8=0;
IF S3BQ1A9A=2 THEN S3BQ1A9A=0;

IF S3BQ1A1=9 THEN S3BQ1A1=.;
IF S3BQ1A2=9 THEN S3BQ1A2=.;
IF S3BQ1A3=9 THEN S3BQ1A3=.;
IF S3BQ1A4=9 THEN S3BQ1A4=.;
IF S3BQ1A5=9 THEN S3BQ1A5=.;
IF S3BQ1A6=9 THEN S3BQ1A6=.;
IF S3BQ1A7=9 THEN S3BQ1A7=.;
IF S3BQ1A8=9 THEN S3BQ1A8=.;
IF S3BQ1A9A=9 THEN S3BQ1A9A=.;

IF S3AQ3B1=9 THEN S3AQ3B1=.;
IF S3AQ3C1=99 THEN S3AQ3C1=.;

IF S2AQ8A=99 THEN S2AQ8A=.;
IF S2AQ8B=99 THEN S2AQ8B=.;

/* DATA BINNING - from QUANTITATIVE to CATEGORICAL BINS */

/*SMOKE = from FREQUENCY IN DAYS converted to NUMBER OF DAYS PER YEAR (days per year) */
IF S3AQ3B1= 1 THEN SMOKEFREQYR = 364;
ELSE IF S3AQ3B1= 2 THEN SMOKEFREQYR = 286;
ELSE IF S3AQ3B1= 3 THEN SMOKEFREQYR = 182;
ELSE IF S3AQ3B1= 4 THEN SMOKEFREQYR = 78;
ELSE IF S3AQ3B1= 5 THEN SMOKEFREQYR = 30;
ELSE IF S3AQ3B1= 6 THEN SMOKEFREQYR = 1;

/* assign NEW VARIABLE pack-years */
CIGSPERYEAR = SMOKEFREQYR * S3AQ3C1;

/* ALCOHOL = from FREQUENCY IN DAYS converted to NUMBER OF DAYS PER YEAR (days per year) */
IF S2AQ8A= 1 THEN ALCOHOLFREQYR = 364;
ELSE IF S2AQ8A= 2 THEN ALCOHOLFREQYR = 286;
ELSE IF S2AQ8A= 3 THEN ALCOHOLFREQYR = 182;
ELSE IF S2AQ8A= 4 THEN ALCOHOLFREQYR = 104;
ELSE IF S2AQ8A= 5 THEN ALCOHOLFREQYR = 52;
ELSE IF S2AQ8A= 6 THEN ALCOHOLFREQYR = 30;
ELSE IF S2AQ8A= 7 THEN ALCOHOLFREQYR = 12;
ELSE IF S2AQ8A= 8 THEN ALCOHOLFREQYR = 9;
ELSE IF S2AQ8A= 9 THEN ALCOHOLFREQYR = 4.5;
ELSE IF S2AQ8A= 10 THEN ALCOHOLFREQYR = 1.5;

/* assign NEW VARIABLE drink-years */
SWIGSPERYEAR = ALCOHOLFREQYR * S2AQ8B;

/*make bins for quantitative explanatory variables cigsperyear */
if CIGSPERYEAR LE 2002 then G4CIGS =1001.5;
else if CIGSPERYEAR LE 4550 then G4CIGS = 3276; 
else if CIGSPERYEAR LE 7644 then G4CIGS =6097; 
else if CIGSPERYEAR GT 7644 then G4CIGS =21658; 

/*make bins for quantitative explanatory variables swigsperyear */

if SWIGSPERYEAR LE 10.5 then G4SWIGS = 6;
else if SWIGSPERYEAR LE 63 then G4SWIGS = 36.75; 
else if SWIGSPERYEAR LE 360 then G4SWIGS =211.5; 
else if SWIGSPERYEAR GT 360 then G4SWIGS =35852; 

/*sort tables by row identifier*/
PROC SORT; by IDNUM;

/*PROC ANOVA; CLASS MAJORDEPLIFE;
MODEL AGE=MAJORDEPLIFE;
MEANS MAJORDEPLIFE /TUKEY;*/

/*PROC ANOVA; CLASS G4CIGS;
MODEL AGE=G4CIGS;
MEANS G4CIGS /TUKEY;*/

/*PROC ANOVA; CLASS G4SWIGS;
MODEL AGE=G4SWIGS;
MEANS G4SWIGS /DUNCAN;*/

/*print each observation with the following variables in columns */
/*PROC PRINT; VAR MAJORDEPLIFE SEX  G4CIGS G4SWIGS S3BQ1A1 S3BQ1A2 S3BQ1A3 S3BQ1A4 S3BQ1A5 S3BQ1A6 S3BQ1A7 S3BQ1A8 S3BQ1A9A*/ ;

/*check center, spread, shape*/
/*PROC UNIVARIATE; VAR  G4CIGS G4SWIGS MAJORDEPLIFE SEX S3BQ1A1 S3BQ1A2 S3BQ1A3 S3BQ1A4 S3BQ1A5 S3BQ1A6 S3BQ1A7 S3BQ1A8 S3BQ1A9A*/ ;

/*check frequency tables*/
/*PROC FREQ; TABLES MAJORDEPLIFE SEX G4CIGS G4SWIGS S2AQ8B S3AQ3C1 S3BQ1A1 S3BQ1A2 S3BQ1A3 S3BQ1A4 S3BQ1A5 S3BQ1A6 S3BQ1A7 S3BQ1A8 S3BQ1A9A*/ ;

RUN;

/*get multiple pairwise tables for chi-squared test*/

/*this is for comparing G4CIGS classes versus majordepressive disorder*/
/*Bonferonni correction is 0.05 divided by 6 equals 0.008333333 */

DATA CIGSCOMPARISON1; set NEW;
IF G4CIGS=1001.5 or G4CIGS=3276;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4CIGS /CHISQ;

RUN;

DATA CIGSCOMPARISON2; set NEW;
IF G4CIGS=1001.5 or G4CIGS=6097;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4CIGS /CHISQ;

RUN;

DATA CIGSCOMPARISON3; set NEW;
IF G4CIGS=1001.5 or G4CIGS=21658;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4CIGS /CHISQ;

RUN;

DATA CIGSCOMPARISON4; set NEW;
IF G4CIGS=3276 or G4CIGS=6097;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4CIGS /CHISQ;

RUN;

DATA CIGSCOMPARISON5; set NEW;
IF G4CIGS=3276 or G4CIGS=21658;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4CIGS /CHISQ;

RUN;

DATA CIGSCOMPARISON6; set NEW;
IF G4CIGS=6097 or G4CIGS=21658;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4CIGS /CHISQ;

RUN;

/*this is for comparing G4SWIGS classes versus majordepressive disorder*/
/*Bonferonni correction is 0.05 divided by 6 equals 0.008333333 */

DATA SWIGSCOMPARISON1; set NEW;
IF G4SWIGS=6 or G4SWIGS=36.75;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4SWIGS /CHISQ;

RUN;

DATA SWIGSCOMPARISON2; set NEW;
IF G4SWIGS=6 or G4SWIGS=211.5;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4SWIGS /CHISQ;

RUN;

DATA SWIGSCOMPARISON3; set NEW;
IF G4SWIGS=6 or G4SWIGS=35852;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4SWIGS /CHISQ;

RUN;

DATA SWIGSCOMPARISON4; set NEW;
IF G4SWIGS=36.75 or G4SWIGS=211.5;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4SWIGS /CHISQ;

RUN;

DATA SWIGSCOMPARISON5; set NEW;
IF G4SWIGS=36.75 or G4SWIGS=35852;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4SWIGS /CHISQ;

RUN;

DATA SWIGSCOMPARISON6; set NEW;
IF G4SWIGS=211.5 or G4SWIGS=35852;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*G4SWIGS /CHISQ;

RUN;

/*COMPARE SUBSTANCES versus MDD */
/*for a 2 by 2 case, we are inclined to call the chi square statistic large, if it is larger than 3.841*/
/* you may check out http://math.hws.edu/javamath/ryan/ChiSquare.html for the table of df versus alpha*/
/*p value is simply 0.05, since it's not a multiple pairwise test*/

DATA COMPARISONS; set NEW;
PROC SORT; by IDNUM;

PROC FREQ; TABLES  MAJORDEPLIFE*S3BQ1A1 /CHISQ; /*sedatives*/
PROC FREQ; TABLES  MAJORDEPLIFE*S3BQ1A2 /CHISQ; /*tranq*/
PROC FREQ; TABLES  MAJORDEPLIFE*S3BQ1A3 /CHISQ; /*opioids*/
PROC FREQ; TABLES  MAJORDEPLIFE*S3BQ1A4 /CHISQ; /*amphet*/
PROC FREQ; TABLES  MAJORDEPLIFE*S3BQ1A5 /CHISQ; /*cannabis*/
PROC FREQ; TABLES  MAJORDEPLIFE*S3BQ1A6 /CHISQ; /*cocaine*/
PROC FREQ; TABLES  MAJORDEPLIFE*S3BQ1A7 /CHISQ; /*hallucinogens*/
PROC FREQ; TABLES  MAJORDEPLIFE*S3BQ1A8 /CHISQ; /*inhalants*/
PROC FREQ; TABLES  MAJORDEPLIFE*S3BQ1A9A /CHISQ; /*heroin*/

RUN;

i have to post this script on this Google blog, because Tumblr does not support the display of colored codes. by the way, this is for my SAS / R course for Wesleyan, and i am cross-referencing this post to my course blog substancestats.tumblr.com.

getwd()

setwd("C:/Users/aseandi/My Library/COURSERA/Passion-Driven Statistics/datasets")

getwd()

# fileurlcsv <- "http://spark-public.s3.amazonaws.com/pdstatistics/data_sets/nesarc_pds.csv"
# download.file(fileurlcsv, destfile="./nesarc_pds.csv")
# list.files(".")
# dateDownloaded <-date()
# dateDownloaded

nesarc <- read.csv("./nesarc_pds.csv")
str(nesarc)
# summary(nesarc)
# sapply(nesarc[1, ], class)
# sum(is.na(nesarc))
# table(is.na(nesarc))

library(Hmisc)
library(RColorBrewer)

library(reshape2)
library(stringr)
library(plyr)

# # --------------------------------------------------------------------------------------------------------
# # # DEFINING THE VARIABLES UNDER STUDY # # #

# study variables: Sex, average daily quantity of alcohol consumed and cigarettes smoked in past 12 months, and use of sedatives, tranquilizers, cannabis, opioids, amphetamines, cocaine, heroine, hallucinogens, and inhalants
# codes: Sex, s2aq8b, s3aq3c1, S3bq1a1, S3bq1a2, S3bq1a3, S3bq1a4, S3bq1a5, S3bq1a6, S3bq1a7, S3bq1a8, S3bq1a9a

#number of alcohol consumed in past 12 months
table(is.na(nesarc$S2AQ8B))
table(nesarc$S2AQ8B)
sum(table(nesarc$S2AQ8B))

#usual quantity of cigarettes smoked
table(is.na(nesarc$S3AQ3C1))
table(nesarc$S3AQ3C1)
sum(table(nesarc$S3AQ3C1))

# # --------------------------------------------------------------------------------------------------------
# # # CREATE NEW DATAFRAME WITH FEWER VARIABLES # # #
mddold <- data.frame(nesarc$MAJORDEPLIFE, nesarc$SEX, nesarc$S3AQ3B1, nesarc$S3AQ3C1, nesarc$S2AQ8A, nesarc$S2AQ8B, nesarc$S3BQ1A1, nesarc$S3BQ1A2, nesarc$S3BQ1A3, nesarc$S3BQ1A4, nesarc$S3BQ1A5, nesarc$S3BQ1A6, nesarc$S3BQ1A7, nesarc$S3BQ1A8, nesarc$S3BQ1A9A)
str(mddold)
names(mddold) <- c("mdd", "sex", "smokefreq", "smoke", "drinkfreq", "alcohol", "sedatives", "tranquilizers", "cannabis", "opioids", "amphetamines", "cocaine", "heroine", "hallucinogens", "inhalants")
str(mddold)

# # --------------------------------------------------------------------------------------------------------
# # # DATA MUNGING # # #
# for the  following variables: "sex", "sedatives", "tranquilizers", "cannabis", "opioids", "amphetamines", "cocaine", "heroine", "hallucinogens", "inhalants"
	# change the binomial setting [1,2] --> [1,0]
	# change [9] to [NA]

# # # DATA MUNGING # # #
# for the  following variables: "alcohol", "smoke", "smokefreq", "drinkfreq"
	# change [99] for smoke quantity, drinks consumed to [NA]
	# change [99] for drinkfreq to [NA]
	# change [9] for smokefreq to [NA]

mddold$sex[mddold$sex==2]=0 # sex=0 female; sex=1 male

mddold$sedatives[mddold$sedatives==2]=0
mddold$tranquilizers[mddold$tranquilizers==2]=0
mddold$cannabis[mddold$cannabis==2]=0
mddold$opioids[mddold$opioids==2]=0
mddold$amphetamines[mddold$amphetamines==2]=0
mddold$cocaine[mddold$cocaine==2]=0
mddold$heroine[mddold$heroine==2]=0
mddold$hallucinogens[mddold$hallucinogens==2]=0
mddold$inhalants[mddold$inhalants==2]=0

mddold$sedatives[mddold$sedatives==9]=NA
mddold$tranquilizers[mddold$tranquilizers==9]=NA
mddold$cannabis[mddold$cannabis==9]=NA
mddold$opioids[mddold$opioids==9]=NA
mddold$amphetamines[mddold$amphetamines==9]=NA
mddold$cocaine[mddold$cocaine==9]=NA
mddold$heroine[mddold$heroine==9]=NA
mddold$hallucinogens[mddold$hallucinogens==9]=NA
mddold$inhalants[mddold$inhalants==9]=NA


mddold$alcohol[mddold$alcohol==99]=NA
mddold$smoke[mddold$smoke==99]=NA
mddold$smokefreq[mddold$smokefreq==9]=NA
mddold$drinkfreq[mddold$drinkfreq==99]=NA

# # # ADD NEW VARIABLES: PACKYEARS and DRINKYEARS # # #

mddold$smokefreqyr[mddold$smokefreq==1]= 364
mddold$smokefreqyr[mddold$smokefreq==2]= 286
mddold$smokefreqyr[mddold$smokefreq==3]= 182
mddold$smokefreqyr[mddold$smokefreq==4]= 78
mddold$smokefreqyr[mddold$smokefreq==5]= 30
mddold$smokefreqyr[mddold$smokefreq==6]= 1
mddold$smokefreqyr[mddold$smokefreq==NA]= NA

mddold$alcoholfreqyr[mddold$drinkfreq==1]=364
mddold$alcoholfreqyr[mddold$drinkfreq==2]=286
mddold$alcoholfreqyr[mddold$drinkfreq==3]=182
mddold$alcoholfreqyr[mddold$drinkfreq==4]=104
mddold$alcoholfreqyr[mddold$drinkfreq==5]=52
mddold$alcoholfreqyr[mddold$drinkfreq==6]=30
mddold$alcoholfreqyr[mddold$drinkfreq==7]=12
mddold$alcoholfreqyr[mddold$drinkfreq==8]=9
mddold$alcoholfreqyr[mddold$drinkfreq==9]=4.5
mddold$alcoholfreqyr[mddold$drinkfreq==10]=1.5
mddold$alcoholfreqyr[mddold$drinkfreq==NA]=NA

mddold$cigsperyear <- (mddold$smokefreqyr * mddold$smoke)
mddold$swigsperyear <- (mddold$alcoholfreqyr * mddold$alcohol)

summary(mddold$smokefreqyr)
summary(mddold$alcoholfreqyr)
summary(mddold$cigsperyear)
summary(mddold$swigsperyear)

table(mddold$smokefreqyr)
table(mddold$alcoholfreqyr)
table(mddold$cigsperyear)
table(mddold$swigsperyear)

mddnew <- mddold
str(mddnew)

# # --------------------------------------------------------------------------------------------------------
# FREQUENCY TABLES
# ("mdd", "sex", "drinkyears", "packyears", "sedatives", "tranquilizers", "cannabis", "opioids", "amphetamines", "cocaine", "heroine", "hallucinogens", "inhalants")

library(gmodels)
CrossTable(mddnew$mdd, mddnew$sex)
CrossTable(mddnew$mdd, mddnew$sedatives)
CrossTable(mddnew$mdd, mddnew$tranquilizers)
CrossTable(mddnew$mdd, mddnew$cannabis)
CrossTable(mddnew$mdd, mddnew$opioids)
CrossTable(mddnew$mdd, mddnew$amphetamines)
CrossTable(mddnew$mdd, mddnew$cocaine)
CrossTable(mddnew$mdd, mddnew$heroine)
CrossTable(mddnew$mdd, mddnew$hallucinogens)
CrossTable(mddnew$mdd, mddnew$inhalants)

library(Hmisc)
mddnew$g4cigs <- cut2(mddnew$cigsperyear, g = 4)
CrossTable(mddnew$mdd, mddnew$g4cigs)
# intervals
# [   1, 2002) | [2002, 4550) | [4550, 7644) | [7644,35672]
  
mddnew$g4swigs <- cut2(mddnew$swigsperyear, g = 4)
CrossTable(mddnew$mdd, mddnew$g4swigs)
# intervals
# [  1.5,   10.5) | [ 10.5,   63.0) | [ 63.0,  360.0) | [360.0,35672.0] 

# # --------------------------------------------------------------------------------------------------------
# # # MULTIVARIATE GRAPHS FOR EXPLORATORY ANALYSIS # # #

library(RColorBrewer)

mypar <- function(a = 1, b = 1, brewer.n = 4, brewer.name = "RdYlGn", ...) {
    par(mar = c(2.5, 2.5, 1.6, 1.1), mgp = c(1.5, 0.5, 0))
    par(mfrow = c(a, b), ...)
    palette(brewer.pal(brewer.n, brewer.name))
}

# create table for mdd vs cigs category
mddvg4cigs = table(mddnew$mdd,mddnew$g4cigs)

# To get the graph we want, we need to exchange the rows in this table
mddvg4cigs = rbind(mddvg4cigs[2,],mddvg4cigs[1,])

# and turn them into percents (dividing by the num. of observations # in each cigs category)
mddvg4cigs[1,]=mddvg4cigs[1,]/table(mddnew$g4cigs)
mddvg4cigs[2,]=mddvg4cigs[2,]/table(mddnew$g4cigs)
str(mddvg4cigs)

# create table for mdd vs alcohol category
mddvg4swigs = table(mddnew$mdd,mddnew$g4swigs)

# To get the graph we want, we need to exchange the rows in this table
mddvg4swigs = rbind(mddvg4swigs[2,],mddvg4swigs[1,])

# and turn them into percents (dividing by the num. of observations # in each swigs category)
mddvg4swigs[1,]=mddvg4swigs[1,]/table(mddnew$g4swigs)
mddvg4swigs[2,]=mddvg4swigs[2,]/table(mddnew$g4swigs)
str(mddvg4swigs)

# create table for mdd vs sex
mddvsex = table(mddnew$mdd,mddnew$sex)

# To get the graph we want, we need to exchange the rows in this table
mddvsex = rbind(mddvsex[2,],mddvsex[1,])

# and turn them into percents (dividing by the num. of observations # in either con or sans sex)
mddvsex[1,]=mddvsex[1,]/table(mddnew$sex)
mddvsex[2,]=mddvsex[2,]/table(mddnew$sex)
str(mddvsex)

# create table for mdd vs sedatives
mddvsedatives = table(mddnew$mdd,mddnew$sedatives)

# To get the graph we want, we need to exchange the rows in this table
mddvsedatives = rbind(mddvsedatives[2,],mddvsedatives[1,])

# and turn them into percents (dividing by the num. of observations # in either con or sans sedatives)
mddvsedatives[1,]=mddvsedatives[1,]/table(mddnew$sedatives)
mddvsedatives[2,]=mddvsedatives[2,]/table(mddnew$sedatives)
str(mddvsedatives)

# create table for mdd vs cannabis
mddvcannabis = table(mddnew$mdd,mddnew$cannabis)

# To get the graph we want, we need to exchange the rows in this table
mddvcannabis = rbind(mddvcannabis[2,],mddvcannabis[1,])

# and turn them into percents (dividing by the num. of observations # in either con or sans cannabis)
mddvcannabis[1,]=mddvcannabis[1,]/table(mddnew$cannabis)
mddvcannabis[2,]=mddvcannabis[2,]/table(mddnew$cannabis)
str(mddvcannabis)

# create table for mdd vs opioids
mddvopioids = table(mddnew$mdd,mddnew$opioids)

# To get the graph we want, we need to exchange the rows in this table
mddvopioids = rbind(mddvopioids[2,],mddvopioids[1,])

# and turn them into percents (dividing by the num. of observations # in either con or sans opioids)
mddvopioids[1,]=mddvopioids[1,]/table(mddnew$opioids)
mddvopioids[2,]=mddvopioids[2,]/table(mddnew$opioids)
str(mddvopioids)

# create table for mdd vs tranquilizers
mddvtranquilizers = table(mddnew$mdd,mddnew$tranquilizers)

# To get the graph we want, we need to exchange the rows in this table
mddvtranquilizers = rbind(mddvtranquilizers[2,],mddvtranquilizers[1,])

# and turn them into percents (dividing by the num. of observations # in either con or sans tranquilizers)
mddvtranquilizers[1,]=mddvtranquilizers[1,]/table(mddnew$tranquilizers)
mddvtranquilizers[2,]=mddvtranquilizers[2,]/table(mddnew$tranquilizers)
str(mddvtranquilizers)

# create table for mdd vs amphetamines
mddvamphetamines = table(mddnew$mdd,mddnew$amphetamines)

# To get the graph we want, we need to exchange the rows in this table
mddvamphetamines = rbind(mddvamphetamines[2,],mddvamphetamines[1,])

# and turn them into percents (dividing by the num. of observations # in either con or sans amphetamines)
mddvamphetamines[1,]=mddvamphetamines[1,]/table(mddnew$amphetamines)
mddvamphetamines[2,]=mddvamphetamines[2,]/table(mddnew$amphetamines)
str(mddvamphetamines)

# create table for mdd vs cocaine
mddvcocaine = table(mddnew$mdd,mddnew$cocaine)

# To get the graph we want, we need to exchange the rows in this table
mddvcocaine = rbind(mddvcocaine[2,],mddvcocaine[1,])

# and turn them into percents (dividing by the num. of observations # in either con or sans cocaine)
mddvcocaine[1,]=mddvcocaine[1,]/table(mddnew$cocaine)
mddvcocaine[2,]=mddvcocaine[2,]/table(mddnew$cocaine)
str(mddvcocaine)

# create table for mdd vs heroine
mddvheroine = table(mddnew$mdd,mddnew$heroine)

# To get the graph we want, we need to exchange the rows in this table
mddvheroine = rbind(mddvheroine[2,],mddvheroine[1,])

# and turn them into percents (dividing by the num. of observations # in either con or sans heroine)
mddvheroine[1,]=mddvheroine[1,]/table(mddnew$heroine)
mddvheroine[2,]=mddvheroine[2,]/table(mddnew$heroine)
str(mddvheroine)

# create table for mdd vs hallucinogens
mddvhallucinogens = table(mddnew$mdd,mddnew$hallucinogens)

# To get the graph we want, we need to exchange the rows in this table
mddvhallucinogens = rbind(mddvhallucinogens[2,],mddvhallucinogens[1,])

# and turn them into percents (dividing by the num. of observations # in either con or sans hallucinogens)
mddvhallucinogens[1,]=mddvhallucinogens[1,]/table(mddnew$hallucinogens)
mddvhallucinogens[2,]=mddvhallucinogens[2,]/table(mddnew$hallucinogens)
str(mddvhallucinogens)

# create table for mdd vs inhalants
mddvinhalants = table(mddnew$mdd,mddnew$inhalants)

# To get the graph we want, we need to exchange the rows in this table
mddvinhalants = rbind(mddvinhalants[2,],mddvinhalants[1,])

# and turn them into percents (dividing by the num. of observations # in either con or sans inhalants)
mddvinhalants[1,]=mddvinhalants[1,]/table(mddnew$inhalants)
mddvinhalants[2,]=mddvinhalants[2,]/table(mddnew$inhalants)
str(mddvinhalants)

mypar(mfrow = c(2,2))
# MDD diagnosis {RESPONSE} by Estimated Cigarette Use per Year {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvg4cigs <- barplot(mddvg4cigs[1,], col=unique(mddnew$g4cigs), xlab="cigarettes per year", ylab="diagnosed depression", cex.axis=0.8)
# MDD diagnosis {RESPONSE} by Estimated Cigarette Use per Year {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvg4cigs <- barplot(mddvg4cigs, col=unique(mddnew$g4cigs), xlab="cigarettes per year", ylab="diagnosed depression", cex.axis=0.8)
# MDD diagnosis {RESPONSE} by Estimated Alcohol Consumed per Year {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvg4swigs <- barplot(mddvg4swigs[1,], col=unique(mddnew$g4swigs), xlab="alcohol per year", ylab="diagnosed depression", cex.axis=0.8)
# MDD diagnosis {RESPONSE} by Estimated Alcohol Consumed per Year {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvg4swigs <- barplot(mddvg4swigs, col=unique(mddnew$g4swigs), xlab="alcohol per year", ylab="diagnosed depression", cex.axis=0.8)
# dev.copy2pdf(file="mdd_cigarettes_alcohol.pdf", height =8, width = 11)

mypar(mfrow = c(2, 5))
# MDD diagnosis {RESPONSE} by Biological Sex {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvsex <- barplot(mddvsex[1,], col=unique(mddnew$sex), xlab="biological sex, 0 - female, 1 - male", ylab="diagnosed depression")
# MDD diagnosis {RESPONSE} by Sedatives Use {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvsedatives <- barplot(mddvsedatives[1,], col=unique(mddnew$sedatives), xlab="sedative use", ylab="diagnosed depression")
# MDD diagnosis {RESPONSE} by Cannabis Use {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvcannabis <- barplot(mddvcannabis[1,], col=unique(mddnew$cannabis), xlab="cannabis use", ylab="diagnosed depression")
# MDD diagnosis {RESPONSE} by tranquilizers use {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvtranquilizers <- barplot(mddvtranquilizers[1,], col=unique(mddnew$tranquilizers), xlab="tranquilizers use", ylab="diagnosed depression")
# MDD diagnosis {RESPONSE} by opioids use {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvopioids <- barplot(mddvopioids[1,], col=unique(mddnew$opioids), xlab="opioids use", ylab="diagnosed depression")
# MDD diagnosis {RESPONSE} by cocaine use {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvcocaine <- barplot(mddvcocaine[1,], col=unique(mddnew$cocaine), xlab="cocaine use", ylab="diagnosed depression")
# MDD diagnosis {RESPONSE} by amphetamines use {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvamphetamines <- barplot(mddvamphetamines[1,], col=unique(mddnew$amphetamines), xlab="amphetamines use", ylab="diagnosed depression")
# MDD diagnosis {RESPONSE} by heroine use {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvheroine <- barplot(mddvheroine[1,], col=unique(mddnew$heroine), xlab="heroine use", ylab="diagnosed depression")
# MDD diagnosis {RESPONSE} by hallucinogens use {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvhallucinogens <- barplot(mddvhallucinogens[1,], col=unique(mddnew$hallucinogens), xlab="hallucinogens use", ylab="diagnosed depression")
# MDD diagnosis {RESPONSE} by inhalants use {EXPLANATORY} among all Adults in the NESARC Study
bp_mddvinhalants <- barplot(mddvinhalants[1,], col=unique(mddnew$inhalants), xlab="inhalants use", ylab="diagnosed depression")
# dev.copy2pdf(file="mdd_sex_substances.pdf", height =8, width = 11)

in any case, i think Google will eventually be able to crawl this R post, and hopefully other beginner programmers like me will find the code useful for their basic graphing.

bovine economics

the world order as explained by two cows

A Tale of Two Cows by kylatomdesign.

becoming the next pope

the short how-to guide for obtaining the title of pope awesome the first

How To Become Pope via Visual.ly

a disney short film

this was the finest 6 minutes in recent animation history. it preceded wreck-it-ralph in most theaters, and the emotions in this short are at par with the full feature film.

best android games of 2012

Here are my personal favorites and what I think are the best Android titles (games category) that were released in 2012.

I would definitely recommend Jetpack Joyride, Clay Jam, Hill Climb Racing, and Zombie Dash in a heartbeat.

Jetpack Joyride	96
Zombie Dash	96
Clay Jam	94
Hill Climb Racing	94
Triple Town	92
Rope Escape	90
Sector Strike	88
iSlash	86
Subway Surfers	86
Music Hero	84
Strikefleet Omega	82
Kick the Boss 2	80
Amazing Alex	75-79
Angry Birds Rio	75-79
Angry Birds Seasons	75-79
Angry Birds Space	75-79
Angry Birds Star Wars	75-79
Bad Piggies	75-79
Chasing Yellow	75-79
Cut the Rope Free	75-79
Icon Pop Quiz	75-79
Radiant Defense	75-79
Samurai vs Zombies	75-79
Swing Shot	75-79
Whale Trail Free	75-79
100 Doors 2013	70-74
Androidify	70-74
Blosics Free	70-74
Flow Free	70-74
Flow Free: Bridges	70-74
GeoQuiz	70-74
GeoQuiz by BrianCafe	70-74
Logo Quiz Philippines	70-74
Math Maniac	70-74
Tetris	70-74
Death Moto	65-69
Drag Racing	65-69
Dragon Fly! Free	65-69
Falling Fred	65-69
Fruit Ninja Free	65-69
Guns N Glory WW2	65-69
Raging Thunder	65-69
Slice It!	65-69
Temple Run	65-69
3D Bowling	60-64
Basketball Mania	60-64
Bubble Blast 2	60-64
Bubble Blast Halloween	60-64
Bubble Worlds	60-64
Burger	60-64
Chess Free	60-64
Jewels Star	60-64
Logo Quiz	60-64
Marble Blast 3	60-64
Paper Toss	60-64
Pinball Pro	60-64
Pool Master Pro	60-64
Random Mahjong	60-64
Unblock Me Free	60-64
Arel Wars2	55-59
Blast Monkeys	55-59
Empire Defense	55-59
Fruit Shoot	55-59
Ninjump	55-59
Rocket Island	55-59
Tank Hero	55-59

Note that most, if not all, of the games listed above have already received excellent user ratings of at least 4 stars out of 5 on Google's Play Store. Personally, I find Google's user scoring system to be a pretty reliable benchmark for Android apps.

In here, I prefer to rank them in terms of replay value, which is really dependent on the challenges and missions that will make you return daily so as to finish the game.

For other players, it could be about the interesting story arc, solid graphics, or the unobtrusive sound that spells the difference.

Some games got lower marks because they constantly need an internet connection, such as ArelWars2.

Other games were also impossibly difficult to finish if you don't prefer purchasing premium items and in-game goods using real money, such as those offered by Samurai vs. Zombies and Radiant Defense. These apps had very stunning graphics and wonderful gameplay, but I was frustated with the constant string of losses due to lack of 'freemium weapons'.

Here, of course, are screenshots of installed games last year.

Since then, I've added quite a few good ones this 2013, most notably a Japanese title called Battle Cats. I think it's going to be a huge hit this year.

Also worth the look is the Android version of Temple Run 2, which already has an Apple version just last week. We'll see how this year unfolds for amateur & professional Android developers.

Navigation

SAS code

R: graphing plots

bovine economics

becoming the next pope

a disney short film

best android games of 2012