###=================================================================### ###=================================================================### ### OMDS Workshop: Exploring Your Data...correlograms ### ### by: Sarah Schwartz ### date: 2-18-16 ###=================================================================### ###=================================================================### ### ----- clear environment ----- ### rm(list = ls()) ### ------------- load packages you expect to use --------------- ### library(haven) # for reading data in from SPSS, SAS, Stata library(ggplot2) # for plotting library(RColorBrewer) # nice looking color palettes library(psych) # lots of cool & helpful functions # For data wrangling library(plyr) # split, apply, combine, arraange, mutate library(dplyr) # subset & join library(tidyr) # gather & spread library(reshape2) # melt & cast ### ------------------ set working directory -----------------------### # this should be a FOLDER, not an specific file (unique for each project) fp = file.path("C:/Users/A00315273/Box Sync", "Office of Research Services", "OMDS", "OMDS Workshops", "OMDS ggplot2") # set the wokring directory, where to get and save 'stuff' setwd(fp) # ask to "get" or show you the working directory getwd() ### ------------------ read in the dataset -----------------------### # the "read_sav" function in the "Haven" package # this reads the SPSS dataset into R's global environment dat = read_sav("IM PA Diet for spaghetti plot.sav") str(dat) # structure, compact display dim(dat) # dimentions, # rows x # columns names(dat) # names of the variables head(dat) # the first few rows (subjects or cases) # frequency of the "Index1" variable table(dat$Index1) # create a new variable called "days" # if Index1 == 1, the days are all zero (its baseline) # if Index1 == 2, the days are v1 to v2 # if Index1 == 3, the days are v1 to v3 dat$days = ifelse(dat$Index1 == 1, 0, ifelse(dat$Index1 == 2, round(dat$daysv1v2), round(dat$daysv1v3))) # check out the distribution in each group describeBy(dat$days, group = dat$Index1) ### ----------------- do some data cleaning ---------------------### # SurveyDate is missing dat[is.na(dat$SurveyDate), c("GMID", "Index1", "days", "daysv1v2", "daysv1v3", "Intrinsic","ModPAminperwk", "SurveyDate")] # remove all the missed assessments dat = dat[!is.na(dat$SurveyDate),] # double check the days variable summary(dat$days) summary(dat$SurveyDate) # change the 'class' of the date variable so it's easy to work with class(dat$SurveyDate) dat$SurveyDate = as.Date(dat$SurveyDate) summary(dat$SurveyDate[dat$Index1 == 1]) # look at all the lines with missing "days" variables dat[is.na(dat$days), c("GMID", "Index1", "days", "daysv1v2", "daysv1v3", "Intrinsic","ModPAminperwk", "SurveyDate")] # check out a specific subject, ID #9059 dat[dat$GMID == 9059, c("GMID", "Index1", "days", "daysv1v2", "daysv1v3", "Intrinsic", "ModPAminperwk", "SurveyDate")] # who is missing "Intrinsic", but has "Moderate activity" dat[is.na(dat$Intrinsic) & !is.na(dat$ModPAminperwk), c("GMID", "Index1", "days", "daysv1v2", "daysv1v3", "Intrinsic", "ModPAminperwk", "SurveyDate")] # who is missing "Moderate activity", but has "Intrinsic" dat[!is.na(dat$Intrinsic) & is.na(dat$ModPAminperwk), c("GMID", "Index1", "days", "daysv1v2", "daysv1v3", "Intrinsic", "ModPAminperwk", "SurveyDate")] # look at all the lines with missing "days" variables dat[is.na(dat$days) & !is.na(dat$SurveyDate), c("GMID", "Index1", "days", "daysv1v2", "daysv1v3", "Intrinsic", "ModPAminperwk", "SurveyDate")] # check out a specific subject, ID #9004 (missed the 1st assessment only) dat[dat$GMID == 9004, c("GMID", "Index1", "days", "daysv1v2", "daysv1v3", "Intrinsic", "ModPAminperwk", "SurveyDate")] # there are some missing "days" variables...not working yet... dat$daysImp = as.numeric(difftime(dat$SurveyDate, as.Date("2014-04-17"), units = "days")) # look at the Imputed days vs. the date difference days dat[1:12, c("GMID", "Index1", "days", "daysImp", "daysv1v2", "daysv1v3", "Intrinsic", "ModPAminperwk", "SurveyDate")] # just look at the ones who were missing their baseline assesment dat[is.na(dat$days) & !is.na(dat$SurveyDate), c("GMID", "Index1", "days", "daysImp", "daysv1v2", "daysv1v3", "Intrinsic", "ModPAminperwk", "SurveyDate")] dat[dat$GMID == 9004, c("GMID", "Index1", "days", "daysImp", "daysv1v2", "daysv1v3", "Intrinsic", "ModPAminperwk", "SurveyDate")] # Use the imputed "days" variable for subjects who missed their baseline assessment dat[is.na(dat$days) & !is.na(dat$SurveyDate), ]$days = dat[is.na(dat$days) & !is.na(dat$SurveyDate), ]$daysImp summary(dat[dat$Index1 == 1, "SurveyDate"], na.rm = TRUE) summary(dat$days) # get ride of times with missing value (moderate physical activity) dat = dat[!is.na(dat$ModPAminperwk), ] ### --------------------- declare factors -----------------------### # set the variables are factors with set labels dat$Index1 = factor(dat$Index1) dat$v1Gender = factor(dat$v1Gender, labels = c("male", "female")) dat$treatmentgrp = factor(dat$treatmentgrp, labels = c("control", "treatment")) #Add labels to a time variable instead of Index number dat$time = factor(dat$Index1, labels = c("Baseline", "Midway", "Finish")) ### ------------------------ pasta time! -------------------------### names(dat) dat[, c("GMID", "treatmentgrp", "v1Gender", "time", "days", "Intrinsic", "ModPAminperwk")] summary(dat$Intrinsic) dat$Intrinsic3 = factor(ntile(dat$Intrinsic, 3)) summary(dat$Intrinsic4) # a first go at the plot - SINGLE PANE SPAGEHTII PLOT ggplot(dat, aes(x = days, y = ModPAminperwk, group = GMID, col = v1Gender)) + geom_point() + geom_line(alpha = 0.2) # a better plot 0 GRID OF SPAGEHTII PLOTS ggplot(filter(dat, !is.na(Intrinsic3)), aes(x = days, y = ModPAminperwk, group = GMID)) + geom_point() + geom_line(alpha = 0.5) + facet_grid(v1Gender ~ Intrinsic3) + coord_cartesian(ylim = c(0, 1000)) # zoom in to see better ggplot(filter(dat, !is.na(Intrinsic3)), aes(x = days, y = ModPAminperwk, group = GMID)) + geom_point() + geom_line(alpha = 0.5) + facet_grid(. ~ Intrinsic3) + coord_cartesian(ylim = c(0, 1000)) ### ------ compare baseline ------ ### ggplot(dat[dat$Index1 == 1, ], aes(x = ModPAminperwk)) + geom_histogram(binwidth = 30) + facet_grid(treatmentgrp ~ .) ggplot(dat[dat$ModPAminperwk <= 500,], aes(x = Index1, y = ModPAminperwk, fill = treatmentgrp)) + geom_boxplot() + theme_bw() # center the time vairable so weeks and weeks-squared aren't correlated dat$cwk = (dat$days - 113)/7 dat$cwkSq = dat$cwk ^2 summary(dat[, c("cwk", "cwkSq")]) # reduce a dataset to include only complete data (mod PA & Intrinsic) dat.full = dat[complete.cases(dat[, c("cwk", "treatmentgrp", "ModPAminperwk", "Intrinsic")]), ] ggplot(dat.full[dat.full$ModPAminperwk <= 1000,], aes(x = treatmentgrp, y = ModPAminperwk, fill = treatmentgrp)) + geom_boxplot() + facet_grid(. ~ time) + theme_bw() # is there a pattern with intrinsic? ggplot(dat.full, aes(x = cwk, y = Intrinsic, group = GMID)) + geom_point() + geom_line() + facet_grid(v1Gender ~ treatmentgrp) + theme_bw()