Reading in and Writing Out Data
#First we need to make a directory for the final project. We can check if it doesn't exist and then make one if it doesn't.
directory_path <- file.path(getwd(), "Final Project")
if (!dir.exists(directory_path)) {
dir.create(directory_path, recursive = TRUE)
}
#Then we can check if our dataset "allkliff" is in the folder and put it in our final project folder if it isn't.
if (!file.exists(file.path(directory_path, "allkliff.sas7bdat"))) {
file.copy("allkliff.sas7bdat", file.path(directory_path, "allkliff.sas7bdat"))
}
setwd(directory_path) #and then set our working directory to this new file
data <- haven::read_sas("allkliff.sas7bdat") #now we can read in the data
str(data) #taking a look at how big the data is: tibble [184 × 19,551] (S3: tbl_df/tbl/data.frame)
## tibble [184 × 19,551] (S3: tbl_df/tbl/data.frame)
## $ FAMID : num [1:184] 500 501 502 503 504 505 506 507 508 509 ...
## ..- attr(*, "label")= chr "FAMID"
## $ age1y : num [1:184] 12.9 13.6 13.4 13.8 13.1 ...
## ..- attr(*, "label")= chr "teens age from access W1 - years"
## $ age2y : num [1:184] 14.3 14.9 15 15.1 14.1 ...
## ..- attr(*, "label")= chr "teens age from access W2 - years"
## $ age3y : num [1:184] 15 16 15.5 15.9 15.1 ...
## ..- attr(*, "label")= chr "teens age from access W3 - years"
## $ age4y : num [1:184] 16.9 17.3 17 17.1 16.4 ...
## ..- attr(*, "label")= chr "teens age from access W4 - years"
## $ age5y : num [1:184] 17.6 18 18 18.1 17.6 ...
## ..- attr(*, "label")= chr "teens age from access W5 - years"
## $ age6y : num [1:184] 19 19 NA 19.5 18.6 ...
## ..- attr(*, "label")= chr "teens age from access W6 - years"
## $ age7y : num [1:184] 19.8 20.4 NA 20.8 19.7 ...
## ..- attr(*, "label")= chr "teens age from access W7 - years"
## $ age8y : num [1:184] 21 NA 21.8 22 21.1 ...
## ..- attr(*, "label")= chr "teens age from access W8 - years"
## $ age9y : num [1:184] 22.6 NA 22.7 22.7 22.5 ...
## ..- attr(*, "label")= chr "teens age from access W9 - years"
## $ age10y : num [1:184] NA 24 23.7 24 23.5 ...
## ..- attr(*, "label")= chr "teens age from access W10 - years"
## $ age11y : num [1:184] 24.7 25.2 24.3 24.7 NA ...
## ..- attr(*, "label")= chr "teens age from access w11 - years"
## $ age12y : num [1:184] 25.4 26.5 25 25.4 25.7 ...
## ..- attr(*, "label")= chr "teens age from access w12 - years"
## $ age13y : num [1:184] 25.9 27.3 26.2 27.4 28.9 ...
## ..- attr(*, "label")= chr "teens age from access w13 - years"
## $ age14y : num [1:184] 26.5 28.1 27 28.2 NA ...
## ..- attr(*, "label")= chr "teens age from access w14 - years"
## $ age15y : num [1:184] 27.6 29.1 28.5 28.8 NA ...
## ..- attr(*, "label")= chr "teens age from access w15 - years"
## $ age16y : num [1:184] 28.5 30.1 NA 29.4 29.4 ...
## ..- attr(*, "label")= chr "teens age from access w16 - years"
## $ age17y : num [1:184] 30 31 NA 31.3 NA ...
## ..- attr(*, "label")= chr "teens age from access w17 - years"
## $ BIRTHDAY : Date[1:184], format: "1985-11-15" "1985-03-19" ...
## $ AGE567y : num [1:184] 21 24 21.8 22 21.1 ...
## ..- attr(*, "label")= chr "teens avg age W5-7"
## $ age18y : num [1:184] 31.6 32.2 NA 32.5 33.1 ...
## ..- attr(*, "label")= chr "teens age from access w18 - years"
## $ age8910y : num [1:184] 21 24 21.8 22 21.1 ...
## ..- attr(*, "label")= chr "teens avg age W9-10"
## $ age1123y : num [1:184] 24.7 25.2 24.3 24.7 25.7 ...
## ..- attr(*, "label")= chr "teens avg age W11-13"
## $ age1456y : num [1:184] 27.6 28.1 NA NA NA ...
## ..- attr(*, "label")= chr "teens avg age W14-16"
## $ age1718y : num [1:184] 30 31 NA NA NA ...
## ..- attr(*, "label")= chr "teens avg age W17-18"
## $ age19y : num [1:184] 33.7 33.9 NA 34.2 34.9 ...
## ..- attr(*, "label")= chr "teens age from access w19 - years"
## $ age20y : num [1:184] NA 35.4 NA 35.5 36 ...
## ..- attr(*, "label")= chr "teens age from access w20 - years"
## $ agey567m : num [1:184] 17.6 18 NA 18.1 17.6 ...
## ..- attr(*, "label")= chr "W5-7 TN Age at Mom Visit"
## $ agey567d : num [1:184] NA NA NA 18.1 17.6 ...
## ..- attr(*, "label")= chr "W5-7 TN Age at Dad Visit"
## $ age21y : num [1:184] NA 36.9 NA 37.1 37.2 ...
## ..- attr(*, "label")= chr "teens age from access w21 - years"
## $ age1920y : num [1:184] NA 33.9 NA 35.5 36 ...
## ..- attr(*, "label")= chr "teens avg age W19-20"
## $ age5r : num [1:184] 17 NA NA NA 16 18 NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age w5 -years"
## $ age6r : num [1:184] NA 16 NA 17 17 26 NA 19 NA 18 ...
## ..- attr(*, "label")= chr "rp age w6 -years"
## $ age7r : num [1:184] 21 NA NA 19 NA 26 NA NA 22 24 ...
## ..- attr(*, "label")= chr "rp age w7 -years"
## $ age8r : num [1:184] 19 NA NA 18 21 27 21 NA 27 25 ...
## ..- attr(*, "label")= chr "rp age w8 -years"
## $ age9r : num [1:184] 21.1 NA NA 22.6 NA ...
## ..- attr(*, "label")= chr "rp age at w9 packet -years"
## $ age10r : num [1:184] NA 20 NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w10 packet -years"
## $ age11r : num [1:184] NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w11 packet -years"
## $ age12r : num [1:184] 23.7 21.4 NA 24.8 25.2 ...
## ..- attr(*, "label")= chr "rp age at w12 packet -years"
## $ age13r : num [1:184] NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w13 packet -years"
## $ age14r : num [1:184] NA 24.1 NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w14 packet -years"
## $ age15r : num [1:184] 22.9 NA NA 28.7 NA ...
## ..- attr(*, "label")= chr "rp age at w15 packet -years"
## $ age16r : num [1:184] NA NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w16 packet -years"
## $ age17r : num [1:184] 25.1 26.8 NA 30 NA ...
## ..- attr(*, "label")= chr "rp age at w17 packet -years"
## $ age18r : num [1:184] NA NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w18 packet -years"
## $ age19r : num [1:184] NA 30.3 NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w19 packet -years"
## $ age20r : num [1:184] NA NA NA 35.3 NA ...
## ..- attr(*, "label")= chr "rp age at w20 packet -years"
## $ age21r : num [1:184] NA NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w21 packet - years"
## $ age8910r : num [1:184] NA 19.9 NA 18.4 21.2 ...
## ..- attr(*, "label")= chr "Age8-10 RP Age"
## $ age1123r : num [1:184] 23.4 20.8 NA 24 25 ...
## ..- attr(*, "label")= chr "Age11-13 RP Age"
## $ age1456r : num [1:184] 22.9 24.1 NA 28.7 NA ...
## $ age1718r : num [1:184] 25.1 26.8 NA 30 NA ...
## $ age1920r : num [1:184] NA 30.3 NA 35.3 NA ...
## $ relstatus5y : num [1:184] 1 NA NA 1 1 1 NA NA NA NA ...
## ..- attr(*, "label")= chr "Rel Status W5: Teen Relationship Status"
## $ relstatus6y : num [1:184] 0 1 NA 1 1 1 0 0 NA 1 ...
## ..- attr(*, "label")= chr "Rel Status W6: Teen Relationship Status"
## $ relstatus7y : num [1:184] 1 0 NA 1 0 1 0 0 1 1 ...
## ..- attr(*, "label")= chr "Rel Status W7: Teen Relationship Status"
## $ relstatus8y : num [1:184] 1 NA NA 1 1 1 1 0 1 1 ...
## ..- attr(*, "label")= chr "Rel Status W8: Teen Relationship Status"
## $ relstatus9y : num [1:184] 1 NA 0 1 1 1 1 0 0 1 ...
## ..- attr(*, "label")= chr "Rel Status W9: Teen Relationship Status"
## $ relstatus10y : num [1:184] NA 1 0 1 0 1 0 1 0 1 ...
## ..- attr(*, "label")= chr "Rel Status W10: Teen Relationship Status"
## $ relstatus11y : num [1:184] 1 1 0 1 NA 1 0 0 0 1 ...
## ..- attr(*, "label")= chr "Rel Status W11: Teen Relationship Status"
## $ relstatus12y : num [1:184] 1 1 0 1 1 1 0 0 0 1 ...
## ..- attr(*, "label")= chr "Rel Status W12: Teen Relationship Status"
## $ relstatus13y : num [1:184] 1 1 0 1 1 1 0 0 0 1 ...
## ..- attr(*, "label")= chr "Rel Status W13: Teen Relationship Status"
## $ relstatus14y : num [1:184] 1 1 0 1 NA 1 0 0 1 1 ...
## ..- attr(*, "label")= chr "Rel Status W14: Teen Relationship Status"
## $ relstatus15y : num [1:184] 1 1 0 1 NA 1 0 0 1 1 ...
## ..- attr(*, "label")= chr "Rel Status W15: Teen Relationship Status"
## $ relstatus16y : num [1:184] 1 1 NA 1 1 1 0 1 1 0 ...
## ..- attr(*, "label")= chr "Rel Status W16: Teen Relationship Status"
## $ relstatus17y : num [1:184] 1 1 NA 1 NA 1 0 1 1 1 ...
## ..- attr(*, "label")= chr "Rel Status W17: Teen Relationship Status"
## $ agefirstsex : num [1:184] 17 13 NA 16 17 16 16 NA 16 15 ...
## ..- attr(*, "label")= chr "teens age when first had sex"
## $ genderid1920r : num [1:184] NA 2 NA 2 NA NA NA 2 1 NA ...
## ..- attr(*, "label")= chr "tf1920r: I identify my gender as"
## $ GND567R : num [1:184] 2 2 NA 2 NA 1 NA NA 1 2 ...
## ..- attr(*, "label")= chr "rp gender w5-7"
## $ GND8910R : num [1:184] 2 2 NA 2 2 1 1 2 NA 2 ...
## ..- attr(*, "label")= chr "rp gender w8-10"
## $ GND1123R : num [1:184] 2 2 NA 2 2 1 NA NA NA 2 ...
## ..- attr(*, "label")= chr "rp gender w11-13"
## $ GND1456R : num [1:184] 2 2 NA 2 NA 1 NA 2 1 2 ...
## ..- attr(*, "label")= chr "rp gender w14-16"
## $ GND1718R : num [1:184] 2 2 NA 2 NA 1 NA 2 NA 2 ...
## ..- attr(*, "label")= chr "rp gender w17-18"
## $ samesex567y : num [1:184] 0 0 NA 0 NA 0 NA NA 0 0 ...
## ..- attr(*, "label")= chr "W567 rp gender dyad, 0= heterosexual, 1= same sex couple"
## $ samesex8910y : num [1:184] 0 0 NA 0 0 0 0 0 NA 0 ...
## ..- attr(*, "label")= chr "W8910 rp gender dyad, 0= heterosexual, 1= same sex couple"
## $ samesex1123y : num [1:184] 0 0 NA 0 0 0 NA NA NA 0 ...
## ..- attr(*, "label")= chr "W1123 rp gender dyad, 0= heterosexual, 1= same sex couple"
## $ samesex1456y : num [1:184] 0 0 NA 0 NA 0 NA 0 0 0 ...
## ..- attr(*, "label")= chr "W1456 rp gender dyad, 0= heterosexual, 1= same sex couple"
## $ samesex1718y : num [1:184] 0 0 NA 0 NA 0 NA 0 NA 0 ...
## ..- attr(*, "label")= chr "W1718 rp gender dyad, 0= heterosexual, 1= same sex couple"
## $ samesex1920y : num [1:184] NA 0 NA 0 NA NA NA 0 0 NA ...
## ..- attr(*, "label")= chr "W1920 rp gender dyad, 0= heterosexual, 1= same sex couple"
## $ oppsexfriend1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W1 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend2y : num [1:184] 0 0 NA 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W2 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend3y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W3 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend4y : num [1:184] 0 0 NA 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W4 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend5y : num [1:184] 0 0 0 0 0 NA 0 0 0 0 ...
## ..- attr(*, "label")= chr "W5 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend6y : num [1:184] 0 0 NA 0 NA NA 0 NA NA 0 ...
## ..- attr(*, "label")= chr "W6 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend7y : num [1:184] 0 0 NA 0 0 NA 1 NA NA NA ...
## ..- attr(*, "label")= chr "W7 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend10y : num [1:184] NA NA 0 0 0 NA 1 0 0 0 ...
## ..- attr(*, "label")= chr "W10 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend11y : num [1:184] 0 0 0 0 NA 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W11 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend12y : num [1:184] 0 0 0 0 NA 0 0 0 0 1 ...
## ..- attr(*, "label")= chr "W12 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend13y : num [1:184] 0 0 NA NA NA 0 NA 0 NA 0 ...
## ..- attr(*, "label")= chr "W13 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend14y : num [1:184] 0 0 0 1 NA 0 1 0 0 0 ...
## ..- attr(*, "label")= chr "W14 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend15y : num [1:184] 0 0 0 1 NA 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W15 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend16y : num [1:184] 0 NA NA 0 0 0 0 0 NA 0 ...
## ..- attr(*, "label")= chr "W16 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend17y : num [1:184] 0 0 NA 0 NA NA 0 0 0 0 ...
## ..- attr(*, "label")= chr "W17 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend18y : num [1:184] 0 0 NA 1 0 0 0 0 0 1 ...
## ..- attr(*, "label")= chr "W18 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend19y : num [1:184] 0 0 NA 1 1 NA 0 0 0 1 ...
## ..- attr(*, "label")= chr "W19 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend20y : num [1:184] NA NA NA 1 NA NA NA NA NA 1 ...
## ..- attr(*, "label")= chr "W20 friendship gender, 1= opposite sex, 0= same sex"
## $ INDID1Y : num [1:184] 1083 1320 1044 1255 1258 ...
## ..- attr(*, "label")= chr "INDID1Y"
## $ GENDER1Y : num [1:184] 1 1 1 1 1 2 2 1 2 1 ...
## ..- attr(*, "label")= chr "DEMO1y: Tn gndr"
## [list output truncated]
head(data) #and then can look at an example of what the data looks like
## # A tibble: 6 × 19,551
## FAMID age1y age2y age3y age4y age5y age6y age7y age8y age9y age10y age11y
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 500 12.9 14.3 15.0 16.9 17.6 19.0 19.8 21.0 22.6 NA 24.7
## 2 501 13.6 14.9 16.0 17.3 18.0 19.0 20.4 NA NA 24.0 25.2
## 3 502 13.4 15 15.5 17.0 18.0 NA NA 21.8 22.7 23.7 24.3
## 4 503 13.8 15.1 15.9 17.1 18.1 19.5 20.8 22.0 22.7 24.0 24.7
## 5 504 13.1 14.1 15.1 16.4 17.6 18.6 19.7 21.1 22.5 23.5 NA
## 6 505 13.0 14.2 NA 16.3 18.1 20.8 20.8 21.4 22.6 23.8 24.3
## # ℹ 19,539 more variables: age12y <dbl>, age13y <dbl>, age14y <dbl>,
## # age15y <dbl>, age16y <dbl>, age17y <dbl>, BIRTHDAY <date>, AGE567y <dbl>,
## # age18y <dbl>, age8910y <dbl>, age1123y <dbl>, age1456y <dbl>,
## # age1718y <dbl>, age19y <dbl>, age20y <dbl>, agey567m <dbl>, agey567d <dbl>,
## # age21y <dbl>, age1920y <dbl>, age5r <dbl>, age6r <dbl>, age7r <dbl>,
## # age8r <dbl>, age9r <dbl>, age10r <dbl>, age11r <dbl>, age12r <dbl>,
## # age13r <dbl>, age14r <dbl>, age15r <dbl>, age16r <dbl>, age17r <dbl>, …
Cleaning Data
# data is very large so we need to take a lot of steps to clean it
selected_columns <- data[, c("FAMID", grep("^[a-zA-Z]+\\d{1,2}[a-zA-Z]$", colnames(data), value = TRUE))] #First thing we can do is take all the variables with the pattern: [letters][1-2 numbers][letters]. For this dataset this represents the variables involving our main participants (rather than their parents, friends, partners, etc).
str(selected_columns) #a little cleaner
## tibble [184 × 15,237] (S3: tbl_df/tbl/data.frame)
## $ FAMID : num [1:184] 500 501 502 503 504 505 506 507 508 509 ...
## ..- attr(*, "label")= chr "FAMID"
## $ age1y : num [1:184] 12.9 13.6 13.4 13.8 13.1 ...
## ..- attr(*, "label")= chr "teens age from access W1 - years"
## $ age2y : num [1:184] 14.3 14.9 15 15.1 14.1 ...
## ..- attr(*, "label")= chr "teens age from access W2 - years"
## $ age3y : num [1:184] 15 16 15.5 15.9 15.1 ...
## ..- attr(*, "label")= chr "teens age from access W3 - years"
## $ age4y : num [1:184] 16.9 17.3 17 17.1 16.4 ...
## ..- attr(*, "label")= chr "teens age from access W4 - years"
## $ age5y : num [1:184] 17.6 18 18 18.1 17.6 ...
## ..- attr(*, "label")= chr "teens age from access W5 - years"
## $ age6y : num [1:184] 19 19 NA 19.5 18.6 ...
## ..- attr(*, "label")= chr "teens age from access W6 - years"
## $ age7y : num [1:184] 19.8 20.4 NA 20.8 19.7 ...
## ..- attr(*, "label")= chr "teens age from access W7 - years"
## $ age8y : num [1:184] 21 NA 21.8 22 21.1 ...
## ..- attr(*, "label")= chr "teens age from access W8 - years"
## $ age9y : num [1:184] 22.6 NA 22.7 22.7 22.5 ...
## ..- attr(*, "label")= chr "teens age from access W9 - years"
## $ age10y : num [1:184] NA 24 23.7 24 23.5 ...
## ..- attr(*, "label")= chr "teens age from access W10 - years"
## $ age11y : num [1:184] 24.7 25.2 24.3 24.7 NA ...
## ..- attr(*, "label")= chr "teens age from access w11 - years"
## $ age12y : num [1:184] 25.4 26.5 25 25.4 25.7 ...
## ..- attr(*, "label")= chr "teens age from access w12 - years"
## $ age13y : num [1:184] 25.9 27.3 26.2 27.4 28.9 ...
## ..- attr(*, "label")= chr "teens age from access w13 - years"
## $ age14y : num [1:184] 26.5 28.1 27 28.2 NA ...
## ..- attr(*, "label")= chr "teens age from access w14 - years"
## $ age15y : num [1:184] 27.6 29.1 28.5 28.8 NA ...
## ..- attr(*, "label")= chr "teens age from access w15 - years"
## $ age16y : num [1:184] 28.5 30.1 NA 29.4 29.4 ...
## ..- attr(*, "label")= chr "teens age from access w16 - years"
## $ age17y : num [1:184] 30 31 NA 31.3 NA ...
## ..- attr(*, "label")= chr "teens age from access w17 - years"
## $ age18y : num [1:184] 31.6 32.2 NA 32.5 33.1 ...
## ..- attr(*, "label")= chr "teens age from access w18 - years"
## $ age19y : num [1:184] 33.7 33.9 NA 34.2 34.9 ...
## ..- attr(*, "label")= chr "teens age from access w19 - years"
## $ age20y : num [1:184] NA 35.4 NA 35.5 36 ...
## ..- attr(*, "label")= chr "teens age from access w20 - years"
## $ age21y : num [1:184] NA 36.9 NA 37.1 37.2 ...
## ..- attr(*, "label")= chr "teens age from access w21 - years"
## $ age5r : num [1:184] 17 NA NA NA 16 18 NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age w5 -years"
## $ age6r : num [1:184] NA 16 NA 17 17 26 NA 19 NA 18 ...
## ..- attr(*, "label")= chr "rp age w6 -years"
## $ age7r : num [1:184] 21 NA NA 19 NA 26 NA NA 22 24 ...
## ..- attr(*, "label")= chr "rp age w7 -years"
## $ age8r : num [1:184] 19 NA NA 18 21 27 21 NA 27 25 ...
## ..- attr(*, "label")= chr "rp age w8 -years"
## $ age9r : num [1:184] 21.1 NA NA 22.6 NA ...
## ..- attr(*, "label")= chr "rp age at w9 packet -years"
## $ age10r : num [1:184] NA 20 NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w10 packet -years"
## $ age11r : num [1:184] NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w11 packet -years"
## $ age12r : num [1:184] 23.7 21.4 NA 24.8 25.2 ...
## ..- attr(*, "label")= chr "rp age at w12 packet -years"
## $ age13r : num [1:184] NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w13 packet -years"
## $ age14r : num [1:184] NA 24.1 NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w14 packet -years"
## $ age15r : num [1:184] 22.9 NA NA 28.7 NA ...
## ..- attr(*, "label")= chr "rp age at w15 packet -years"
## $ age16r : num [1:184] NA NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w16 packet -years"
## $ age17r : num [1:184] 25.1 26.8 NA 30 NA ...
## ..- attr(*, "label")= chr "rp age at w17 packet -years"
## $ age18r : num [1:184] NA NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w18 packet -years"
## $ age19r : num [1:184] NA 30.3 NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w19 packet -years"
## $ age20r : num [1:184] NA NA NA 35.3 NA ...
## ..- attr(*, "label")= chr "rp age at w20 packet -years"
## $ age21r : num [1:184] NA NA NA NA NA ...
## ..- attr(*, "label")= chr "rp age at w21 packet - years"
## $ relstatus5y : num [1:184] 1 NA NA 1 1 1 NA NA NA NA ...
## ..- attr(*, "label")= chr "Rel Status W5: Teen Relationship Status"
## $ relstatus6y : num [1:184] 0 1 NA 1 1 1 0 0 NA 1 ...
## ..- attr(*, "label")= chr "Rel Status W6: Teen Relationship Status"
## $ relstatus7y : num [1:184] 1 0 NA 1 0 1 0 0 1 1 ...
## ..- attr(*, "label")= chr "Rel Status W7: Teen Relationship Status"
## $ relstatus8y : num [1:184] 1 NA NA 1 1 1 1 0 1 1 ...
## ..- attr(*, "label")= chr "Rel Status W8: Teen Relationship Status"
## $ relstatus9y : num [1:184] 1 NA 0 1 1 1 1 0 0 1 ...
## ..- attr(*, "label")= chr "Rel Status W9: Teen Relationship Status"
## $ relstatus10y : num [1:184] NA 1 0 1 0 1 0 1 0 1 ...
## ..- attr(*, "label")= chr "Rel Status W10: Teen Relationship Status"
## $ relstatus11y : num [1:184] 1 1 0 1 NA 1 0 0 0 1 ...
## ..- attr(*, "label")= chr "Rel Status W11: Teen Relationship Status"
## $ relstatus12y : num [1:184] 1 1 0 1 1 1 0 0 0 1 ...
## ..- attr(*, "label")= chr "Rel Status W12: Teen Relationship Status"
## $ relstatus13y : num [1:184] 1 1 0 1 1 1 0 0 0 1 ...
## ..- attr(*, "label")= chr "Rel Status W13: Teen Relationship Status"
## $ relstatus14y : num [1:184] 1 1 0 1 NA 1 0 0 1 1 ...
## ..- attr(*, "label")= chr "Rel Status W14: Teen Relationship Status"
## $ relstatus15y : num [1:184] 1 1 0 1 NA 1 0 0 1 1 ...
## ..- attr(*, "label")= chr "Rel Status W15: Teen Relationship Status"
## $ relstatus16y : num [1:184] 1 1 NA 1 1 1 0 1 1 0 ...
## ..- attr(*, "label")= chr "Rel Status W16: Teen Relationship Status"
## $ relstatus17y : num [1:184] 1 1 NA 1 NA 1 0 1 1 1 ...
## ..- attr(*, "label")= chr "Rel Status W17: Teen Relationship Status"
## $ oppsexfriend1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W1 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend2y : num [1:184] 0 0 NA 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W2 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend3y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W3 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend4y : num [1:184] 0 0 NA 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W4 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend5y : num [1:184] 0 0 0 0 0 NA 0 0 0 0 ...
## ..- attr(*, "label")= chr "W5 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend6y : num [1:184] 0 0 NA 0 NA NA 0 NA NA 0 ...
## ..- attr(*, "label")= chr "W6 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend7y : num [1:184] 0 0 NA 0 0 NA 1 NA NA NA ...
## ..- attr(*, "label")= chr "W7 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend10y : num [1:184] NA NA 0 0 0 NA 1 0 0 0 ...
## ..- attr(*, "label")= chr "W10 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend11y : num [1:184] 0 0 0 0 NA 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W11 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend12y : num [1:184] 0 0 0 0 NA 0 0 0 0 1 ...
## ..- attr(*, "label")= chr "W12 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend13y : num [1:184] 0 0 NA NA NA 0 NA 0 NA 0 ...
## ..- attr(*, "label")= chr "W13 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend14y : num [1:184] 0 0 0 1 NA 0 1 0 0 0 ...
## ..- attr(*, "label")= chr "W14 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend15y : num [1:184] 0 0 0 1 NA 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W15 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend16y : num [1:184] 0 NA NA 0 0 0 0 0 NA 0 ...
## ..- attr(*, "label")= chr "W16 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend17y : num [1:184] 0 0 NA 0 NA NA 0 0 0 0 ...
## ..- attr(*, "label")= chr "W17 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend18y : num [1:184] 0 0 NA 1 0 0 0 0 0 1 ...
## ..- attr(*, "label")= chr "W18 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend19y : num [1:184] 0 0 NA 1 1 NA 0 0 0 1 ...
## ..- attr(*, "label")= chr "W19 friendship gender, 1= opposite sex, 0= same sex"
## $ oppsexfriend20y : num [1:184] NA NA NA 1 NA NA NA NA NA 1 ...
## ..- attr(*, "label")= chr "W20 friendship gender, 1= opposite sex, 0= same sex"
## $ INDID1Y : num [1:184] 1083 1320 1044 1255 1258 ...
## ..- attr(*, "label")= chr "INDID1Y"
## $ GENDER1Y : num [1:184] 1 1 1 1 1 2 2 1 2 1 ...
## ..- attr(*, "label")= chr "DEMO1y: Tn gndr"
## $ INDID16Y : num [1:184] 1083 1320 NA 1255 1258 ...
## ..- attr(*, "label")= chr "Teen's Indid:"
## $ STATU16Y : num [1:184] 1 1 NA 1 1 1 1 1 1 1 ...
## ..- attr(*, "label")= chr "Teen's Status ID:"
## $ gender15y : num [1:184] 1 1 1 1 1 2 2 1 2 1 ...
## ..- attr(*, "label")= chr "Current Gender W15"
## $ gender16y : num [1:184] 1 1 1 1 1 2 2 1 2 1 ...
## ..- attr(*, "label")= chr "Current Gender W16"
## $ gender17y : num [1:184] 1 1 1 1 1 2 2 1 2 1 ...
## ..- attr(*, "label")= chr "Current Gender W17"
## $ gender18y : num [1:184] 1 1 1 1 1 2 2 1 2 1 ...
## ..- attr(*, "label")= chr "Current Gender W18"
## $ gender19y : num [1:184] 1 1 1 1 6 2 2 1 2 1 ...
## ..- attr(*, "label")= chr "Current Gender W19"
## $ gender20y : num [1:184] 1 1 1 1 5 2 2 1 2 1 ...
## ..- attr(*, "label")= chr "Current Gender W20"
## $ gender21y : num [1:184] 1 1 1 1 5 2 2 1 2 1 ...
## ..- attr(*, "label")= chr "Current Gender W21"
## $ gender22y : num [1:184] 1 1 1 1 5 2 2 1 2 1 ...
## ..- attr(*, "label")= chr "Current Gender W22"
## $ cpage1c : num [1:184] 13.3 14 13.9 13.9 14 ...
## ..- attr(*, "label")= chr "cp age at w1 visit -years"
## $ cpage2c : num [1:184] 15.1 15.6 15.1 15 14.8 ...
## ..- attr(*, "label")= chr "cp age at w2 visit -years"
## $ cpage3c : num [1:184] 15.9 16.5 16.6 16.1 16.2 ...
## ..- attr(*, "label")= chr "cp age at w3 visit -years"
## $ cpage4c : num [1:184] 18.1 17.5 NA 17.2 17.5 ...
## ..- attr(*, "label")= chr "cp age at w4 visit -years"
## $ cpage5c : num [1:184] 19.5 18.2 17.4 18.1 18.4 ...
## ..- attr(*, "label")= chr "cp age at w5 visit -years"
## $ momage1m : num [1:184] 33.6 40.7 44.9 44.9 46.4 ...
## ..- attr(*, "label")= chr "moms age from access W1 - years"
## $ momage4m : num [1:184] 37.6 44.4 48.5 48.1 49.7 ...
## ..- attr(*, "label")= chr "moms age from access W4 - years"
## $ momage5m : num [1:184] 38.2 45.1 49.5 49.1 50.9 ...
## ..- attr(*, "label")= chr "moms age from access W5 - years"
## $ momage6m : num [1:184] 39.6 46.1 NA 50.5 52 ...
## ..- attr(*, "label")= chr "moms age from access W6 - years"
## $ momage7m : num [1:184] 40.5 47.5 NA 51.8 53 ...
## ..- attr(*, "label")= chr "moms age from access W7 - years"
## $ dadage1d : num [1:184] NA 38.2 47.5 49.4 45.6 ...
## ..- attr(*, "label")= chr "dads age from access W1 - years"
## $ dadage4d : num [1:184] NA 41.9 51.2 52.6 48.9 ...
## ..- attr(*, "label")= chr "dads age from access W4 - years"
## $ dadage5d : num [1:184] NA 42.6 52.2 53.6 50.1 ...
## ..- attr(*, "label")= chr "dads age from access W5 - years"
## $ dadage6d : num [1:184] NA 43.6 NA 55 51.2 ...
## ..- attr(*, "label")= chr "dads age from access W6 - years"
## $ dadage7d : num [1:184] NA 45 NA 56.3 52.2 ...
## ..- attr(*, "label")= chr "dads age from access W7 - years"
## $ mfriends1y : num [1:184] 12 NA 12 12 10 2 0 12 0 5 ...
## ..- attr(*, "label")= chr "CC1Y: teen number of male friends"
## $ ffriends1y : num [1:184] 0 NA 0 0 2 10 12 0 12 5 ...
## ..- attr(*, "label")= chr "CC1Y: teen number of female friends"
## [list output truncated]
#need to narrow down the dataset even more so let's only look at variables collected in value one (values <2), but make sure to keep the id column
numeric_values <- sapply(names(selected_columns)[-1], function(x) as.numeric(gsub("\\D", "", x)))
selected_columns <- selected_columns[, c(TRUE, numeric_values < 2)]
#we've narrowed down the dataset from 19k to 1.5k but that's still pretty large. Let's look at variables that are reported by or observed of our orginial participants (variable ends in y)
selected_cols = vector(mode = "character", length = 0)
#loop through column names
for (col in colnames(selected_columns)) {
#check if column name ends with "y"
if (endsWith(col, "y")) {
#add column name to selected_cols vector
selected_cols <- c(selected_cols, col)
}
}
#select the columns we found
selected_columns <- selected_columns[, selected_cols]
#we've narrowed it down to 437. But this is still too large. Since we only have 184 observations, we can make our analysis later easier if we just throw out columns that don't have that much data. Let's set a threshold of about no more than 50% empty
#calculate the proportion of NA values in each column
na_prop <- colMeans(is.na(selected_columns))
#threshold for proportion of NA values
threshold <- 0.5
#keep columns where proportion of NA values is below the threshold
selected_columns <- selected_columns[, na_prop < threshold]
str(selected_columns) #narrowed it down a little more
## tibble [184 × 398] (S3: tbl_df/tbl/data.frame)
## $ age1y : num [1:184] 12.9 13.6 13.4 13.8 13.1 ...
## ..- attr(*, "label")= chr "teens age from access W1 - years"
## $ oppsexfriend1y: num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "W1 friendship gender, 1= opposite sex, 0= same sex"
## $ mfriends1y : num [1:184] 12 NA 12 12 10 2 0 12 0 5 ...
## ..- attr(*, "label")= chr "CC1Y: teen number of male friends"
## $ ffriends1y : num [1:184] 0 NA 0 0 2 10 12 0 12 5 ...
## ..- attr(*, "label")= chr "CC1Y: teen number of female friends"
## $ numfriends1y : num [1:184] 12 NA 12 12 12 12 12 12 12 12 ...
## ..- attr(*, "label")= chr "CC1Y= teen number of friends on CC"
## $ propmf1y : num [1:184] 1 NA 1 1 0.833 ...
## ..- attr(*, "label")= chr "CC1Y= proportion of male friends"
## $ propff1y : num [1:184] 0 NA 0 0 0.167 ...
## ..- attr(*, "label")= chr "CC1Y= proportion of female friends"
## $ oppsexf1y : num [1:184] 0 NA 0 0 2 2 0 0 0 5 ...
## ..- attr(*, "label")= chr "CC1Y= teen number of oppostie sex friends"
## $ oppsexbf1y : num [1:184] 0 NA 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "CC1Y= teen best friend opposite sex, 1 = yes"
## $ propoppsexf1y : num [1:184] 0 NA 0 0 0.167 ...
## ..- attr(*, "label")= chr "CC1Y= proportion of opposite sex friends"
## $ ring1y : num [1:184] 4 NA 4 4 4 4 4 4 4 4 ...
## ..- attr(*, "label")= chr "CC1Y = number of friends in inner circle"
## $ cengc1y : num [1:184] 3.25 3 3.33 2.88 2.88 ...
## ..- attr(*, "label")= chr "SBTY1C CP seeking support: youths engagement"
## $ cgisc1y : num [1:184] 3.25 3.25 3.5 2.88 3.25 ...
## ..- attr(*, "label")= chr "SBTY1C CP seeking support: youth given-instrumental support"
## $ cgesc1y : num [1:184] 0 2.5 1.333 0.125 0.5 ...
## ..- attr(*, "label")= chr "SBTY1C CP seeking support: youth given-emotional support"
## $ cintc1y : num [1:184] 4 4 4 4 4 3.75 4 3.25 3.75 3.5 ...
## ..- attr(*, "label")= chr "SBTY1C CP seeking support: youth interpretation"
## $ emp1y : num [1:184] 2.62 2.88 3 2.25 2.69 ...
## ..- attr(*, "label")= chr "SBTy1c: Jessie Empathy Measure"
## $ perpov1y : num [1:184] 0.152 2.128 3.04 4.255 3.636 ...
## ..- attr(*, "label")= chr "perpov1y: W1 percentage of poverty line based on household size and income"
## $ aaqagg1y : num [1:184] 27 28 32 30 21 43 37 20 29 39 ...
## ..- attr(*, "label")= chr "AAQ1y attitudes toward aggression"
## $ legit1y : num [1:184] 6 8 13 11 8 16 13 8 10 15 ...
## ..- attr(*, "label")= chr "AAQ1y legitamacy of aggression"
## $ esteem1y : num [1:184] 5 7 3 3 3 8 3 3 5 6 ...
## ..- attr(*, "label")= chr "AAQ1y aggression increases slf estm"
## $ negimg1y : num [1:184] 5 5 7 6 4 6 8 3 5 8 ...
## ..- attr(*, "label")= chr "AAQ1y aggress decr neg imge"
## $ deserv1y : num [1:184] 4 3 5 6 3 7 7 3 3 5 ...
## ..- attr(*, "label")= chr "AAQ1y victims dsrv aggression"
## $ suffer1y : num [1:184] 7 5 4 4 3 6 6 3 6 5 ...
## ..- attr(*, "label")= chr "AAQ1y: victims do not suffer"
## $ smokev1y : num [1:184] 0 1 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr smkd cig"
## $ oftcig1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Freq smk pst 30 days"
## $ alcoev1y : num [1:184] 0 1 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr drnk alc"
## $ oftalc1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Freq drnk past 30 days"
## $ drkalc1y : num [1:184] 0 NA 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: How oftn drunk"
## $ marjev1y : num [1:184] 0 0 0 0 0 0 0 0 0 1 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr try marijuana"
## $ oftmar1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: How oftn use marijuana"
## $ hallev1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr try hallucinogens"
## $ barbev1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr try barbituates"
## $ trnkev1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr try tranquilizers"
## $ amphev1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr try amphetamines"
## $ inhaev1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr try inhalents"
## $ heroev1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr try heroine"
## $ cocaev1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr try cocaine"
## $ othrev1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr try other drug"
## $ nmeoth1y : chr [1:184] ".................................................." ".................................................." ".................................................." ".................................................." ...
## ..- attr(*, "label")= chr "ALC1Y: Name of other drug tried"
## $ hdrgev1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ACL1y: Total hard drugs tried"
## $ totsoft1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## $ softev1y : num [1:184] 0 1 0 0 0 0 0 0 0 1 ...
## ..- attr(*, "label")= chr "ALC1Y: Ever used soft drugs"
## $ oftsoft1y : num [1:184] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Freq soft drugs 30 days"
## $ apibase1y : num [1:184] 25 35 32 31 34 31 21 32 37 30 ...
## ..- attr(*, "label")= chr "API1Y: baseline score"
## $ apicomp1y : num [1:184] 38 35 28 37 32 31 40 30 40 39 ...
## ..- attr(*, "label")= chr "API1Y: competent influence items score"
## $ apiinco1y : num [1:184] 35 31 32 33 24 20 18 21 38 8 ...
## ..- attr(*, "label")= chr "API1Y: incompetent influence items score"
## $ apiirr1y : num [1:184] 61 61 56 61 43 47 35 54 64 46 ...
## ..- attr(*, "label")= chr "API1Y: Teen response to peer pressure to do something irresponsible"
## $ apiagg1y : num [1:184] 17 15 18 16 9 13 12 27 5 24 ...
## ..- attr(*, "label")= chr "API1Y: Teen response to peer aggression"
## $ apiinfcmp1y : num [1:184] 13 0 -4 6 -2 0 19 -2 3 9 ...
## ..- attr(*, "label")= chr "API1Y: change from competent influence"
## $ apiinficmp1y : num [1:184] 10 -4 0 2 -10 -11 -3 -11 1 -22 ...
## ..- attr(*, "label")= chr "API1Y: change from incompetent influence"
## $ apimse1y : num [1:184] 3 4.4 4.4 4.4 4.4 4.4 3.6 4.6 4.8 4.2 ...
## ..- attr(*, "label")= chr "APISE1Y: mean self-efficacy score"
## $ rsnsd1y : num [1:184] NA 2 2.38 1 3 ...
## $ confd1y : num [1:184] NA 2.25 2.12 2.5 3.75 ...
## $ colld1y : num [1:184] NA 1 0.625 0 0 ...
## $ blrsd1y : num [1:184] NA 0 0 0 0.25 ...
## $ prssd1y : num [1:184] NA 1.5 0 0 0 0 0 0.5 NA NA ...
## $ querd1y : num [1:184] NA 1.5 2.38 0 0 ...
## $ valdd1y : num [1:184] NA 1 0.75 2 1.25 0.5 1.75 0.5 NA NA ...
## $ enggd1y : num [1:184] NA 2.5 3 2.75 3.5 ...
## $ intrd1y : num [1:184] NA 0.75 0.5 1.5 1 ...
## $ critd1y : num [1:184] NA 0 0 0 0 0 0 0 NA NA ...
## $ psatd1y : num [1:184] NA 2.12 2.25 1.75 3.38 ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Pos Auto to Dad"
## $ psrld1y : num [1:184] NA 1.75 1.88 2.38 2.38 ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Pos Rel to Dad"
## $ ngatd1y : num [1:184] NA 0.8333 0.2083 0 0.0833 ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Neg Auto to Dad"
## $ ngrld1y : num [1:184] NA 0.375 0.25 0.75 0.5 ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Neg Rel to Dad"
## $ psard1y : num [1:184] NA 1.94 2.06 2.06 2.88 ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Pos A/R to Dad"
## $ ngard1y : num [1:184] NA 0.65 0.225 0.3 0.25 0.1 0.15 0.7 NA NA ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Neg A/R to Dad"
## $ rsnsm1y : num [1:184] 1.25 1.25 2 2 4 2.5 3.6 1.5 1.25 2.25 ...
## $ confm1y : num [1:184] 1.5 2 2.5 2.5 4 3.5 3.8 1 1.5 2.5 ...
## $ collm1y : num [1:184] 1 0 0 1 0 0 0.1 0.25 1.75 0.25 ...
## $ blrsm1y : num [1:184] 0.5 0 0 0 0.75 0 2.2 0 1 1 ...
## $ prssm1y : num [1:184] 0.25 0.25 0 1.75 0 0.25 0.7 0 0.5 1.25 ...
## $ querm1y : num [1:184] 0.75 0.75 0 0.75 0 0 1.3 0 1 1.25 ...
## $ valdm1y : num [1:184] 0 1 1 1.75 0.75 0.75 0.8 0.5 0.25 1.25 ...
## $ enggm1y : num [1:184] 1.25 2.5 3 2.5 3 2.5 2.2 1.5 2 2.5 ...
## $ intrm1y : num [1:184] 2.75 0.25 0.75 0.75 0 0.75 1.2 0.75 0.25 1.25 ...
## $ critm1y : num [1:184] 0.75 0 0 0 0 0 0.2 0 0.25 1 ...
## $ psatm1y : num [1:184] 1.38 1.62 2.25 2.25 4 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Pos Auto to Mom"
## $ psrlm1y : num [1:184] 0.625 1.75 2 2.125 1.875 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Pos Rel to Mom"
## $ ngatm1y : num [1:184] 0.5833 0.0833 0 0.9167 0.25 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Neg Auto to Mom"
## $ ngrlm1y : num [1:184] 1.75 0.125 0.375 0.375 0 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Neg Rel to Mom"
## $ psarm1y : num [1:184] 1 1.69 2.12 2.19 2.94 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Pos A/R to Mom"
## $ ngarm1y : num [1:184] 1.05 0.1 0.15 0.7 0.15 0.2 0.88 0.2 0.75 0.95 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Neg A/R to Mom"
## $ sxasbk1y : num [1:184] 9 10 10 10 10 7 8 10 9 8 ...
## ..- attr(*, "label")= chr "ASBK1Y: Attitiudes abt Age Having sex"
## $ romasbk1y : num [1:184] 28 11 17 22 17 19 22 15 18 21 ...
## ..- attr(*, "label")= chr "ASBK1Y: Attitudes abt Age Romantic acts"
## $ sumasbk1y : num [1:184] 37 21 27 32 27 26 30 25 27 29 ...
## ..- attr(*, "label")= chr "ASBK1Y: Sum all attitudes questions"
## $ cdisum1y : num [1:184] 2.86 0 4.38 2.62 4.38 ...
## ..- attr(*, "label")= chr "CDI1Y: summary score"
## $ cesumk1y : num [1:184] 23.1 28 27 30 24 ...
## ..- attr(*, "label")= chr "CESBK1Y: CESBQ Summary Rejection Score"
## $ cesumm1y : num [1:184] 15 26 25 23 28.9 ...
## ..- attr(*, "label")= chr "CESBM1Y: CESBQ Summary Score"
## $ acptd1y : num [1:184] NA 27 27 19 25 30 21 25 NA 15 ...
## ..- attr(*, "label")= chr "CRPD1Y: Acceptance v. Rejection"
## $ pctrld1y : num [1:184] 11.2 17 17.8 13 11 ...
## ..- attr(*, "label")= chr "CRPD1Y: Psych. Control v. Autonomy"
## $ fctrld1y : num [1:184] NA 20 15.6 23 14 ...
## ..- attr(*, "label")= chr "CRPD1Y: Firm v. Lax Control"
## $ acptm1y : num [1:184] 26 29 27 28 28 21 24 23 21 26 ...
## ..- attr(*, "label")= chr "CRPM1Y: Acceptance v. Rejection"
## $ pctrlm1y : num [1:184] 21 19 15 11 16 14 14 12 17 14 ...
## ..- attr(*, "label")= chr "CRPM1Y: Psych. Control v. Autonomy"
## $ fctrlm1y : num [1:184] 21 19 17 20 18 23 17 21 20 19 ...
## ..- attr(*, "label")= chr "CRPM1Y: Firm v. Lax Control"
## $ emotab1y : num [1:184] 6 5 5 6 5 6 7 6 7 9 ...
## ..- attr(*, "label")= chr "CTQ1Y: Emotional abuse scale"
## $ emotne1y : num [1:184] 9 5 5 8 6 6 9 7 12 9 ...
## ..- attr(*, "label")= chr "CTQ1Y: Emotional neglect scale"
## $ minden1y : num [1:184] 14 12 10 10 9 13 9 13 7 11 ...
## ..- attr(*, "label")= chr "CTQ1Y: Minimization denial scale"
## $ ygiveind1y : num [1:184] 2 1 0 1 1 0 3 1 NA 1 ...
## ..- attr(*, "label")= chr "CTS tn rpt dad gives in-W1"
## [list output truncated]
#just to see what we're working with at this point let's look at the counts for na value
colSums(!is.na(selected_columns)) #and they all look like they have a good amount of data
## age1y oppsexfriend1y mfriends1y ffriends1y numfriends1y
## 184 180 183 183 183
## propmf1y propff1y oppsexf1y oppsexbf1y propoppsexf1y
## 183 183 183 183 183
## ring1y cengc1y cgisc1y cgesc1y cintc1y
## 183 172 172 172 172
## emp1y perpov1y aaqagg1y legit1y esteem1y
## 165 183 183 184 183
## negimg1y deserv1y suffer1y smokev1y oftcig1y
## 183 183 184 183 183
## alcoev1y oftalc1y drkalc1y marjev1y oftmar1y
## 183 183 149 183 183
## hallev1y barbev1y trnkev1y amphev1y inhaev1y
## 183 182 183 183 183
## heroev1y cocaev1y othrev1y nmeoth1y hdrgev1y
## 183 183 183 184 183
## totsoft1y softev1y oftsoft1y apibase1y apicomp1y
## 183 184 184 172 172
## apiinco1y apiirr1y apiagg1y apiinfcmp1y apiinficmp1y
## 172 172 172 172 172
## apimse1y rsnsd1y confd1y colld1y blrsd1y
## 148 97 97 97 97
## prssd1y querd1y valdd1y enggd1y intrd1y
## 97 97 97 97 97
## critd1y psatd1y psrld1y ngatd1y ngrld1y
## 97 97 97 97 97
## psard1y ngard1y rsnsm1y confm1y collm1y
## 97 97 166 166 166
## blrsm1y prssm1y querm1y valdm1y enggm1y
## 166 166 166 166 166
## intrm1y critm1y psatm1y psrlm1y ngatm1y
## 166 166 166 166 166
## ngrlm1y psarm1y ngarm1y sxasbk1y romasbk1y
## 166 166 166 180 181
## sumasbk1y cdisum1y cesumk1y cesumm1y acptd1y
## 181 182 183 182 156
## pctrld1y fctrld1y acptm1y pctrlm1y fctrlm1y
## 158 157 182 182 182
## emotab1y emotne1y minden1y ygiveind1y ycomprd1y
## 181 183 181 158 156
## dignored1y yignored1y dcontd1y ycontd1y dhurttotd1y
## 158 155 158 158 150
## dadposd1y dadvrbd1y dadphyd1y dadphwd1y dadphld1y
## 153 158 150 150 150
## dadabusd1y tndposd1y tndvrbd1y tndphyd1y tndphwd1y
## 158 152 158 150 150
## tndphld1y tndabusd1y ygiveinm1y mgiveinm1y ymcomprm1y
## 150 158 176 177 176
## yignorem1y mignorem1y ycontm1y mcontm1y mhurttotm1y
## 176 176 177 177 172
## momposm1y momvrbm1y momphym1y momphwm1y momphlm1y
## 180 182 176 176 176
## momabusm1y tnmposm1y tnmvrbm1y tnmphym1y tnmphwm1y
## 182 180 174 176 176
## tnmphlm1y tnmabusm1y ego1y cmpanc1y resoc1y
## 176 182 182 183 177
## helpc1y xchngc1y btrayc1y validc1y supportc1y
## 183 183 183 182 181
## totfqqc1y romapp1y slfwth1y schcmp1y clsfrd1y
## 182 174 182 182 179
## socacc1y drgalc1y totprs1y intprs1y extprs1y
## 182 178 183 183 181
## doprs1y areprs1y socprs1y lead1y follow1y
## 182 183 181 181 183
## agghart1y formagg1y relagg1y ovrtagg1y soccom1y
## 181 179 181 181 177
## trustd1y commd1y aliend1y tattd1y trustk1y
## 157 158 158 157 183
## commk1y alienk1y tattk1y trustm1y commm1y
## 183 183 183 182 182
## alienm1y tattm1y parnmdis1y Parlosper1y Parlosrper1y
## 181 181 179 179 179
## Parlossper1y eidsum1y eachev1y eblong1y eotror1y
## 179 184 183 184 184
## ebehav1y ethslf1y eslfop1y ethmom1y ethdad1y
## 184 184 184 182 156
## expach1y explrn1y affblg1y affirm1y explor1y
## 183 184 184 184 184
## forcls1y minrty1y srda1y lsrda1y srdb1y
## 184 184 183 183 183
## lsrdb1y srdt1y srdvi1y lsrdvi1y srdnv1y
## 183 183 183 183 183
## lsrdnv1y pqcba1y pqdta1y pqhga1y pqtma1y
## 183 180 182 182 181
## pqada1y pqnfta1y pqnsma1y pqgra1y pqfma1y
## 182 182 180 180 182
## pqwta1y pqntsa1y pqlka1y pqncta1y pqtoa1y
## 164 182 182 182 182
## pqrla1y pqsca1y pqcbk1y pqdtk1y pqhgk1y
## 182 182 181 180 182
## pqtmk1y pqadk1y pqnftk1y pqftk1y pqnsmk1y
## 182 182 182 180 174
## pqgrk1y pqfmk1y pqwtk1y pqtsk1y pqlkk1y
## 181 177 179 175 181
## pqctk1y pqnctk1y pqtok1y pqrlk1y pqsck1y
## 180 179 181 181 181
## yposc1y ynegc1y yengc1y ydomc1y ytlkc1y
## 179 179 179 179 179
## yslfc1y ycisc1y ycesc1y ysatc1y posm1y
## 179 179 179 179 168
## negm1y engm1y domm1y tlkm1y slfm1y
## 168 168 168 168 168
## cism1y cesm1y satm1y mwrm1y mitm1y
## 168 168 168 168 168
## mcom1y gone1y alntgt1y bgfrnd1y bgtotl1y
## 168 183 135 183 171
## gngnow1y lnggng1y hhands1y ksslps1y mkeout1y
## 183 183 182 180 178
## prsrks1y sexevr1y consex1y prtctn1y condom1y
## 183 179 180 184 184
## prsusx1y noprsu1y prsosx1y noprso1y frndsx1y
## 184 164 183 184 182
## frndpr1y tlkphn1y hngsch1y hnghse1y hngout1y
## 174 137 138 138 138
## eatlnc1y schlev1y famact1y hwtgtr1y gotgtr1y
## 138 138 138 138 138
## atlkto1y adoing1y atlkme1y ahvfun1y acares1y
## 184 184 184 184 184
## acareo1y aohvbg1y ahelps1y alearn1y apoplr1y
## 184 184 184 184 184
## dlesfr1y djelos1y dtimet1y dexpns1y dngsme1y
## 184 184 184 184 184
## dstbrn1y dbossy1y dfight1y dworry1y gendis1y
## 184 184 184 184 184
## expdate1y dating1y risksex1y pref1y impact1y
## 183 183 184 184 184
## like1y dlike1y sgroup1y tphgk1y tptrk1y
## 184 184 184 184 184
## tpnruk1y tpndrk1y tpdrk1y tpnsxk1y tpsxk1y
## 183 183 183 182 183
## tptok1y tpgenk1y tpnnbk1y tphga1y tptra1y
## 184 183 182 180 180
## tprua1y tpndra1y tpnsxa1y tptoa1y ysrag1y
## 180 180 179 180 183
## ysrdel1y ysrhos1y ysrhyp1y ysrim1y extrn1y
## 183 183 182 183 183
## tnconf1y miscond1y academ1y popular1y tnconfk1y
## 183 183 181 181 176
## miscondk1y academk1y populark1y wqtoa1y wpadva1y
## 176 173 174 166 172
## wptkba1y wptkga1y wpmoda1y wppuna1y wphga1y
## 173 181 172 164 167
## wptra1y wprua1y wpnsxa1y wpndra1y rsnsc1y
## 171 168 167 170 178
## confc1y avoidc1y ovperc1y prssc1y warmc1y
## 178 178 178 178 178
## collabc1y rudec1y tdom1y psatc1y ngatc1y
## 178 178 178 178 178
## psrlc1y ngrlc1y psarc1y ngarc1y wqtok1y
## 178 178 178 178 118
## wqadvk1y wqtkbk1y wqtkgk1y wqmodk1y wqpunk1y
## 182 182 182 182 182
## wqhgk1y wqtrk1y wqnruk1y wqnsxk1y wqsxk1y
## 183 183 183 141 144
## wqndrk1y wqdrk1y wqgenk1y wqnnbk1y pat1y
## 182 135 120 180 169
## cposc1y cnegc1y ctlkc1y cslfc1y bmi1y
## 172 172 172 172 179
## gpa1y intdate1y intdatec1y
## 145 184 184
#ok let's narrow it down even further (we still have almost 400 variables) by getting rid of variables that are very skewed or have excess kurtosis
#let's keep all the columns that are 'good" (skew between -2 and 2 and kurtosis between -7 and 7) in these vectors:
good_column_names <- c()
good_column_indices <- c()
#loop through each column in the dataframe
for (i in seq_along(selected_columns)) {
#check if the column contains numeric values (they all should)
if (is.numeric(selected_columns[[i]])) {
#remove NA values from the column
cleaned_column <- na.omit(selected_columns[[i]])
#check if the column has non-zero variance (error encountered if we don't do this)
if (var(cleaned_column) > 0) {
#calculate skewness and kurtosis for the cleaned column
skewness_value <- skew(cleaned_column)
kurtosis_value <- kurtosi(cleaned_column)
#check if skewness and kurtosis are within -2 and 2 or -7 and 7
if (!is.nan(skewness_value) && !is.nan(kurtosis_value) &&
abs(skewness_value) <= 2 && abs(kurtosis_value) <= 7) {
#ff yes, store the column name and index
good_column_names <- c(good_column_names, names(selected_columns)[i])
good_column_indices <- c(good_column_indices, i)
} else {
#print information about columns with potential issues (just a check to make sure everything looks good and we can potentially include columns that are close to the cutoffs)
cat(sprintf("Column: %s\nSkewness: %.2f\nKurtosis: %.2f\n\n",
names(selected_columns)[i], skewness_value, kurtosis_value))
}
} else {
#print message if the column has zero variance (just a check)
cat(sprintf("Column: %s has zero variance.\n\n", names(selected_columns)[i]))
}
}
}
## Column: oppsexfriend1y
## Skewness: 13.19
## Kurtosis: 173.03
##
## Column: numfriends1y
## Skewness: -2.90
## Kurtosis: 8.30
##
## Column: oppsexbf1y
## Skewness: 3.02
## Kurtosis: 7.18
##
## Column: ring1y
## Skewness: -5.75
## Kurtosis: 31.25
##
## Column: oftcig1y
## Skewness: 3.89
## Kurtosis: 13.18
##
## Column: oftalc1y
## Skewness: 5.98
## Kurtosis: 41.36
##
## Column: drkalc1y
## Skewness: 6.98
## Kurtosis: 52.85
##
## Column: marjev1y
## Skewness: 3.16
## Kurtosis: 8.03
##
## Column: oftmar1y
## Skewness: 8.13
## Kurtosis: 73.14
##
## Column: hallev1y
## Skewness: 9.33
## Kurtosis: 85.54
##
## Column: barbev1y
## Skewness: 5.73
## Kurtosis: 31.05
##
## Column: trnkev1y
## Skewness: 13.31
## Kurtosis: 176.03
##
## Column: amphev1y
## Skewness: 13.31
## Kurtosis: 176.03
##
## Column: inhaev1y
## Skewness: 7.55
## Kurtosis: 55.37
##
## Column: heroev1y has zero variance.
##
## Column: cocaev1y has zero variance.
##
## Column: othrev1y
## Skewness: 9.33
## Kurtosis: 85.54
##
## Column: hdrgev1y
## Skewness: 6.84
## Kurtosis: 55.12
##
## Column: totsoft1y
## Skewness: 7.84
## Kurtosis: 75.07
##
## Column: oftsoft1y
## Skewness: 4.96
## Kurtosis: 28.48
##
## Column: colld1y
## Skewness: 2.57
## Kurtosis: 7.11
##
## Column: prssd1y
## Skewness: 2.24
## Kurtosis: 6.75
##
## Column: critd1y
## Skewness: 4.49
## Kurtosis: 21.29
##
## Column: critm1y
## Skewness: 3.86
## Kurtosis: 19.45
##
## Column: emotab1y
## Skewness: 2.20
## Kurtosis: 5.23
##
## Column: dhurttotd1y
## Skewness: 2.21
## Kurtosis: 4.78
##
## Column: dadphyd1y
## Skewness: 3.05
## Kurtosis: 8.89
##
## Column: dadphwd1y
## Skewness: 3.20
## Kurtosis: 9.78
##
## Column: dadphld1y
## Skewness: 2.22
## Kurtosis: 3.80
##
## Column: dadabusd1y
## Skewness: 2.92
## Kurtosis: 6.95
##
## Column: tndphyd1y
## Skewness: 2.42
## Kurtosis: 5.52
##
## Column: tndphwd1y
## Skewness: 4.36
## Kurtosis: 24.26
##
## Column: tndabusd1y
## Skewness: 2.93
## Kurtosis: 6.85
##
## Column: mhurttotm1y
## Skewness: 2.27
## Kurtosis: 4.65
##
## Column: momphym1y
## Skewness: 6.94
## Kurtosis: 58.50
##
## Column: momphwm1y
## Skewness: 11.47
## Kurtosis: 139.98
##
## Column: momphlm1y
## Skewness: 2.75
## Kurtosis: 9.84
##
## Column: momabusm1y
## Skewness: 3.20
## Kurtosis: 10.17
##
## Column: tnmphym1y
## Skewness: 4.97
## Kurtosis: 27.79
##
## Column: tnmphwm1y
## Skewness: 7.51
## Kurtosis: 65.57
##
## Column: tnmphlm1y
## Skewness: 3.92
## Kurtosis: 16.49
##
## Column: tnmabusm1y
## Skewness: 5.44
## Kurtosis: 28.14
##
## Column: srda1y
## Skewness: 5.89
## Kurtosis: 47.49
##
## Column: lsrda1y
## Skewness: 3.37
## Kurtosis: 17.78
##
## Column: srdb1y
## Skewness: 5.96
## Kurtosis: 49.36
##
## Column: lsrdb1y
## Skewness: 3.33
## Kurtosis: 17.06
##
## Column: srdt1y
## Skewness: 3.97
## Kurtosis: 20.86
##
## Column: srdvi1y
## Skewness: 4.02
## Kurtosis: 22.85
##
## Column: lsrdvi1y
## Skewness: 2.29
## Kurtosis: 6.95
##
## Column: srdnv1y
## Skewness: 7.12
## Kurtosis: 57.09
##
## Column: lsrdnv1y
## Skewness: 5.20
## Kurtosis: 34.12
##
## Column: pqwta1y
## Skewness: 2.06
## Kurtosis: 3.26
##
## Column: pqftk1y
## Skewness: 3.43
## Kurtosis: 11.69
##
## Column: pqwtk1y
## Skewness: 2.14
## Kurtosis: 3.42
##
## Column: pqctk1y
## Skewness: 3.90
## Kurtosis: 14.26
##
## Column: ynegc1y
## Skewness: 3.36
## Kurtosis: 15.23
##
## Column: negm1y
## Skewness: 2.97
## Kurtosis: 10.39
##
## Column: lnggng1y
## Skewness: 11.72
## Kurtosis: 147.07
##
## Column: sexevr1y
## Skewness: 2.98
## Kurtosis: 6.91
##
## Column: consex1y
## Skewness: 3.13
## Kurtosis: 7.82
##
## Column: prtctn1y
## Skewness: -9.36
## Kurtosis: 86.04
##
## Column: condom1y
## Skewness: -3.32
## Kurtosis: 9.10
##
## Column: noprsu1y
## Skewness: 7.29
## Kurtosis: 62.05
##
## Column: prsosx1y
## Skewness: 7.55
## Kurtosis: 55.37
##
## Column: noprso1y
## Skewness: 13.34
## Kurtosis: 177.03
##
## Column: frndpr1y
## Skewness: 4.85
## Kurtosis: 32.74
##
## Column: tlkphn1y has zero variance.
##
## Column: gotgtr1y
## Skewness: 2.75
## Kurtosis: 5.59
##
## Column: atlkto1y
## Skewness: -4.44
## Kurtosis: 17.82
##
## Column: adoing1y
## Skewness: -2.11
## Kurtosis: 2.46
##
## Column: alearn1y
## Skewness: 2.69
## Kurtosis: 5.24
##
## Column: apoplr1y
## Skewness: 4.79
## Kurtosis: 21.06
##
## Column: dlesfr1y
## Skewness: -2.91
## Kurtosis: 6.49
##
## Column: dngsme1y
## Skewness: 2.49
## Kurtosis: 4.24
##
## Column: dstbrn1y
## Skewness: 3.32
## Kurtosis: 9.10
##
## Column: dbossy1y
## Skewness: 4.79
## Kurtosis: 21.06
##
## Column: dfight1y
## Skewness: 5.77
## Kurtosis: 31.45
##
## Column: dworry1y
## Skewness: 6.51
## Kurtosis: 40.55
##
## Column: gendis1y
## Skewness: 9.36
## Kurtosis: 86.04
##
## Column: risksex1y
## Skewness: 3.58
## Kurtosis: 11.91
##
## Column: dlike1y
## Skewness: 2.55
## Kurtosis: 10.20
##
## Column: tpnruk1y
## Skewness: 2.67
## Kurtosis: 7.28
##
## Column: tpdrk1y
## Skewness: 2.50
## Kurtosis: 6.28
##
## Column: tpsxk1y
## Skewness: 2.02
## Kurtosis: 3.85
##
## Column: ysrhos1y
## Skewness: 2.32
## Kurtosis: 8.26
##
## Column: rudec1y
## Skewness: 2.02
## Kurtosis: 4.43
##
## Column: ngrlc1y
## Skewness: 2.02
## Kurtosis: 4.43
##
## Column: wqnruk1y
## Skewness: 2.10
## Kurtosis: 4.53
##
## Column: wqsxk1y
## Skewness: 2.38
## Kurtosis: 5.67
##
## Column: wqdrk1y
## Skewness: 2.79
## Kurtosis: 7.39
##
## Column: cnegc1y
## Skewness: 2.68
## Kurtosis: 7.07
selected_columns <- selected_columns[, good_column_names] #now that we've figured out the weird columns, let's only keep the ones with good skew and kurtosis
str(selected_columns) #we've narrowed it down to a good amount of variables now
## tibble [184 × 305] (S3: tbl_df/tbl/data.frame)
## $ age1y : num [1:184] 12.9 13.6 13.4 13.8 13.1 ...
## ..- attr(*, "label")= chr "teens age from access W1 - years"
## $ mfriends1y : num [1:184] 12 NA 12 12 10 2 0 12 0 5 ...
## ..- attr(*, "label")= chr "CC1Y: teen number of male friends"
## $ ffriends1y : num [1:184] 0 NA 0 0 2 10 12 0 12 5 ...
## ..- attr(*, "label")= chr "CC1Y: teen number of female friends"
## $ propmf1y : num [1:184] 1 NA 1 1 0.833 ...
## ..- attr(*, "label")= chr "CC1Y= proportion of male friends"
## $ propff1y : num [1:184] 0 NA 0 0 0.167 ...
## ..- attr(*, "label")= chr "CC1Y= proportion of female friends"
## $ oppsexf1y : num [1:184] 0 NA 0 0 2 2 0 0 0 5 ...
## ..- attr(*, "label")= chr "CC1Y= teen number of oppostie sex friends"
## $ propoppsexf1y: num [1:184] 0 NA 0 0 0.167 ...
## ..- attr(*, "label")= chr "CC1Y= proportion of opposite sex friends"
## $ cengc1y : num [1:184] 3.25 3 3.33 2.88 2.88 ...
## ..- attr(*, "label")= chr "SBTY1C CP seeking support: youths engagement"
## $ cgisc1y : num [1:184] 3.25 3.25 3.5 2.88 3.25 ...
## ..- attr(*, "label")= chr "SBTY1C CP seeking support: youth given-instrumental support"
## $ cgesc1y : num [1:184] 0 2.5 1.333 0.125 0.5 ...
## ..- attr(*, "label")= chr "SBTY1C CP seeking support: youth given-emotional support"
## $ cintc1y : num [1:184] 4 4 4 4 4 3.75 4 3.25 3.75 3.5 ...
## ..- attr(*, "label")= chr "SBTY1C CP seeking support: youth interpretation"
## $ emp1y : num [1:184] 2.62 2.88 3 2.25 2.69 ...
## ..- attr(*, "label")= chr "SBTy1c: Jessie Empathy Measure"
## $ perpov1y : num [1:184] 0.152 2.128 3.04 4.255 3.636 ...
## ..- attr(*, "label")= chr "perpov1y: W1 percentage of poverty line based on household size and income"
## $ aaqagg1y : num [1:184] 27 28 32 30 21 43 37 20 29 39 ...
## ..- attr(*, "label")= chr "AAQ1y attitudes toward aggression"
## $ legit1y : num [1:184] 6 8 13 11 8 16 13 8 10 15 ...
## ..- attr(*, "label")= chr "AAQ1y legitamacy of aggression"
## $ esteem1y : num [1:184] 5 7 3 3 3 8 3 3 5 6 ...
## ..- attr(*, "label")= chr "AAQ1y aggression increases slf estm"
## $ negimg1y : num [1:184] 5 5 7 6 4 6 8 3 5 8 ...
## ..- attr(*, "label")= chr "AAQ1y aggress decr neg imge"
## $ deserv1y : num [1:184] 4 3 5 6 3 7 7 3 3 5 ...
## ..- attr(*, "label")= chr "AAQ1y victims dsrv aggression"
## $ suffer1y : num [1:184] 7 5 4 4 3 6 6 3 6 5 ...
## ..- attr(*, "label")= chr "AAQ1y: victims do not suffer"
## $ smokev1y : num [1:184] 0 1 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr smkd cig"
## $ alcoev1y : num [1:184] 0 1 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "label")= chr "ALC1Y: Evr drnk alc"
## $ softev1y : num [1:184] 0 1 0 0 0 0 0 0 0 1 ...
## ..- attr(*, "label")= chr "ALC1Y: Ever used soft drugs"
## $ apibase1y : num [1:184] 25 35 32 31 34 31 21 32 37 30 ...
## ..- attr(*, "label")= chr "API1Y: baseline score"
## $ apicomp1y : num [1:184] 38 35 28 37 32 31 40 30 40 39 ...
## ..- attr(*, "label")= chr "API1Y: competent influence items score"
## $ apiinco1y : num [1:184] 35 31 32 33 24 20 18 21 38 8 ...
## ..- attr(*, "label")= chr "API1Y: incompetent influence items score"
## $ apiirr1y : num [1:184] 61 61 56 61 43 47 35 54 64 46 ...
## ..- attr(*, "label")= chr "API1Y: Teen response to peer pressure to do something irresponsible"
## $ apiagg1y : num [1:184] 17 15 18 16 9 13 12 27 5 24 ...
## ..- attr(*, "label")= chr "API1Y: Teen response to peer aggression"
## $ apiinfcmp1y : num [1:184] 13 0 -4 6 -2 0 19 -2 3 9 ...
## ..- attr(*, "label")= chr "API1Y: change from competent influence"
## $ apiinficmp1y : num [1:184] 10 -4 0 2 -10 -11 -3 -11 1 -22 ...
## ..- attr(*, "label")= chr "API1Y: change from incompetent influence"
## $ apimse1y : num [1:184] 3 4.4 4.4 4.4 4.4 4.4 3.6 4.6 4.8 4.2 ...
## ..- attr(*, "label")= chr "APISE1Y: mean self-efficacy score"
## $ rsnsd1y : num [1:184] NA 2 2.38 1 3 ...
## $ confd1y : num [1:184] NA 2.25 2.12 2.5 3.75 ...
## $ blrsd1y : num [1:184] NA 0 0 0 0.25 ...
## $ querd1y : num [1:184] NA 1.5 2.38 0 0 ...
## $ valdd1y : num [1:184] NA 1 0.75 2 1.25 0.5 1.75 0.5 NA NA ...
## $ enggd1y : num [1:184] NA 2.5 3 2.75 3.5 ...
## $ intrd1y : num [1:184] NA 0.75 0.5 1.5 1 ...
## $ psatd1y : num [1:184] NA 2.12 2.25 1.75 3.38 ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Pos Auto to Dad"
## $ psrld1y : num [1:184] NA 1.75 1.88 2.38 2.38 ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Pos Rel to Dad"
## $ ngatd1y : num [1:184] NA 0.8333 0.2083 0 0.0833 ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Neg Auto to Dad"
## $ ngrld1y : num [1:184] NA 0.375 0.25 0.75 0.5 ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Neg Rel to Dad"
## $ psard1y : num [1:184] NA 1.94 2.06 2.06 2.88 ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Pos A/R to Dad"
## $ ngard1y : num [1:184] NA 0.65 0.225 0.3 0.25 0.1 0.15 0.7 NA NA ...
## ..- attr(*, "label")= chr "W1 Dad AR: Youths Neg A/R to Dad"
## $ rsnsm1y : num [1:184] 1.25 1.25 2 2 4 2.5 3.6 1.5 1.25 2.25 ...
## $ confm1y : num [1:184] 1.5 2 2.5 2.5 4 3.5 3.8 1 1.5 2.5 ...
## $ collm1y : num [1:184] 1 0 0 1 0 0 0.1 0.25 1.75 0.25 ...
## $ blrsm1y : num [1:184] 0.5 0 0 0 0.75 0 2.2 0 1 1 ...
## $ prssm1y : num [1:184] 0.25 0.25 0 1.75 0 0.25 0.7 0 0.5 1.25 ...
## $ querm1y : num [1:184] 0.75 0.75 0 0.75 0 0 1.3 0 1 1.25 ...
## $ valdm1y : num [1:184] 0 1 1 1.75 0.75 0.75 0.8 0.5 0.25 1.25 ...
## $ enggm1y : num [1:184] 1.25 2.5 3 2.5 3 2.5 2.2 1.5 2 2.5 ...
## $ intrm1y : num [1:184] 2.75 0.25 0.75 0.75 0 0.75 1.2 0.75 0.25 1.25 ...
## $ psatm1y : num [1:184] 1.38 1.62 2.25 2.25 4 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Pos Auto to Mom"
## $ psrlm1y : num [1:184] 0.625 1.75 2 2.125 1.875 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Pos Rel to Mom"
## $ ngatm1y : num [1:184] 0.5833 0.0833 0 0.9167 0.25 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Neg Auto to Mom"
## $ ngrlm1y : num [1:184] 1.75 0.125 0.375 0.375 0 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Neg Rel to Mom"
## $ psarm1y : num [1:184] 1 1.69 2.12 2.19 2.94 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Pos A/R to Mom"
## $ ngarm1y : num [1:184] 1.05 0.1 0.15 0.7 0.15 0.2 0.88 0.2 0.75 0.95 ...
## ..- attr(*, "label")= chr "ARM1Y Youths Neg A/R to Mom"
## $ sxasbk1y : num [1:184] 9 10 10 10 10 7 8 10 9 8 ...
## ..- attr(*, "label")= chr "ASBK1Y: Attitiudes abt Age Having sex"
## $ romasbk1y : num [1:184] 28 11 17 22 17 19 22 15 18 21 ...
## ..- attr(*, "label")= chr "ASBK1Y: Attitudes abt Age Romantic acts"
## $ sumasbk1y : num [1:184] 37 21 27 32 27 26 30 25 27 29 ...
## ..- attr(*, "label")= chr "ASBK1Y: Sum all attitudes questions"
## $ cdisum1y : num [1:184] 2.86 0 4.38 2.62 4.38 ...
## ..- attr(*, "label")= chr "CDI1Y: summary score"
## $ cesumk1y : num [1:184] 23.1 28 27 30 24 ...
## ..- attr(*, "label")= chr "CESBK1Y: CESBQ Summary Rejection Score"
## $ cesumm1y : num [1:184] 15 26 25 23 28.9 ...
## ..- attr(*, "label")= chr "CESBM1Y: CESBQ Summary Score"
## $ acptd1y : num [1:184] NA 27 27 19 25 30 21 25 NA 15 ...
## ..- attr(*, "label")= chr "CRPD1Y: Acceptance v. Rejection"
## $ pctrld1y : num [1:184] 11.2 17 17.8 13 11 ...
## ..- attr(*, "label")= chr "CRPD1Y: Psych. Control v. Autonomy"
## $ fctrld1y : num [1:184] NA 20 15.6 23 14 ...
## ..- attr(*, "label")= chr "CRPD1Y: Firm v. Lax Control"
## $ acptm1y : num [1:184] 26 29 27 28 28 21 24 23 21 26 ...
## ..- attr(*, "label")= chr "CRPM1Y: Acceptance v. Rejection"
## $ pctrlm1y : num [1:184] 21 19 15 11 16 14 14 12 17 14 ...
## ..- attr(*, "label")= chr "CRPM1Y: Psych. Control v. Autonomy"
## $ fctrlm1y : num [1:184] 21 19 17 20 18 23 17 21 20 19 ...
## ..- attr(*, "label")= chr "CRPM1Y: Firm v. Lax Control"
## $ emotne1y : num [1:184] 9 5 5 8 6 6 9 7 12 9 ...
## ..- attr(*, "label")= chr "CTQ1Y: Emotional neglect scale"
## $ minden1y : num [1:184] 14 12 10 10 9 13 9 13 7 11 ...
## ..- attr(*, "label")= chr "CTQ1Y: Minimization denial scale"
## $ ygiveind1y : num [1:184] 2 1 0 1 1 0 3 1 NA 1 ...
## ..- attr(*, "label")= chr "CTS tn rpt dad gives in-W1"
## $ ycomprd1y : num [1:184] 3 4 1 2 1 1 1 1 NA 3 ...
## ..- attr(*, "label")= chr "CTS tn rpt dad we compromise-W1"
## $ dignored1y : num [1:184] 1 2 0 1 1 1 2 1 NA 1 ...
## ..- attr(*, "label")= chr "CTS tn rpt dad ignores-W1"
## $ yignored1y : num [1:184] 3 1 0 2 1 0 1 0 NA 0 ...
## ..- attr(*, "label")= chr "CTS tn rpt dad: I ignore-W1"
## $ dcontd1y : num [1:184] 1 1 2 1 2 0 1 3 NA 3 ...
## ..- attr(*, "label")= chr "CTS tn rpt dad: I continue-W1"
## $ ycontd1y : num [1:184] 1 2 0 2 0 0 1 3 NA 3 ...
## ..- attr(*, "label")= chr "CTS tn rpt dad: dad continues-W1"
## $ dadposd1y : num [1:184] 5 6 5 4 8 7 12 9 NA 7 ...
## ..- attr(*, "label")= chr "CTS teen rpt dad->tn pos reasoning-W1"
## $ dadvrbd1y : num [1:184] 3 3 2 5 2 1 1 0 NA 6 ...
## ..- attr(*, "label")= chr "CTS teen rpt dad->tn psych aggr-W1"
## $ tndposd1y : num [1:184] 4 5 8 5 8 8 12 8 NA 8 ...
## ..- attr(*, "label")= chr "CTS teen rpt tn->dad pos reasoning-W1"
## $ tndvrbd1y : num [1:184] 9 9 10 3 2 3 14 7 NA 7 ...
## ..- attr(*, "label")= chr "CTS teen rpt tn->dad psych aggr-W1"
## $ tndphld1y : num [1:184] NA NA NA NA NA NA 0 NA NA NA ...
## ..- attr(*, "label")= chr "CTS teen rpt tn->dad WGT LOG phy aggr-W1"
## $ ygiveinm1y : num [1:184] 2 2 2 1 1 1 1 1 2 2 ...
## ..- attr(*, "label")= chr "CTS tn rpt mom - i give in - w1"
## $ mgiveinm1y : num [1:184] 2 2 0 2 1 1 3 3 2 2 ...
## ..- attr(*, "label")= chr "CTS tn rpt mom - she gives in - w1"
## $ ymcomprm1y : num [1:184] 3 4 4 3 1 0 2 3 2 3 ...
## ..- attr(*, "label")= chr "CTS tn rpt mom - we compromise - w1"
## $ yignorem1y : num [1:184] 1 1 1 1 1 1 2 1 2 1 ...
## ..- attr(*, "label")= chr "CTS tn rpt mom - i ignore - w1"
## $ mignorem1y : num [1:184] 1 0 0 1 1 0 0 3 1 1 ...
## ..- attr(*, "label")= chr "CTS tn report mom - she ignores - w1"
## $ ycontm1y : num [1:184] 1 0 2 1 1 1 3 1 1 1 ...
## ..- attr(*, "label")= chr "CTS tn rpt mom - i continue - w1"
## $ mcontm1y : num [1:184] 1 0 0 0 3 0 0 1 1 1 ...
## ..- attr(*, "label")= chr "CTS tn rpt mom - she continues - w1"
## $ momposm1y : num [1:184] 2 4 4 6 4 0 10 6 6 8 ...
## ..- attr(*, "label")= chr "CTS teen rpt mom->tn pos reasoning-W1"
## $ momvrbm1y : num [1:184] 5 2 1.2 0 9 1 1 0 5 1 ...
## ..- attr(*, "label")= chr "CTS teen rpt mom->tn psych aggr-W1"
## $ tnmposm1y : num [1:184] 4 4 8 5 8 0 9 6 7 8 ...
## ..- attr(*, "label")= chr "CTS teen rpt tn->mom pos reasoning-W1"
## $ tnmvrbm1y : num [1:184] 1 2 11 2 4 1 10 1 7 1 ...
## ..- attr(*, "label")= chr "CTS teen rpt tn->mom psych aggr-W1"
## $ ego1y : num [1:184] 3.9 4.02 4.48 4.5 4.65 ...
## $ cmpanc1y : num [1:184] 16 15 22 23 19 25 22 23 16 24 ...
## ..- attr(*, "label")= chr "FQQC1Y: Companionship and Recreation"
## $ resoc1y : num [1:184] 5 12 14 15 15 14 13 11 6 10 ...
## ..- attr(*, "label")= chr "FQQC1Y: Conflict Resolution"
## $ helpc1y : num [1:184] 26 31 38 31 27 41 36 36 30 31 ...
## ..- attr(*, "label")= chr "FQQC1Y: Help and Guidance"
## $ xchngc1y : num [1:184] 10 13 21 17 19 25 20 19 21 15 ...
## ..- attr(*, "label")= chr "FQQC1Y: Intimate Exchange"
## [list output truncated]
df = selected_columns #so let's save that into a separate dataframe
write.csv(df, "new_data.csv", row.names = FALSE) #let's save this data into a new csv
#lets make a function to save the summary stats we want to look at later
#define the function to calculate summary statistics
summary_stats <- function(df) {
#extract numeric columns (excluding the first column)
numeric_cols <- df[, sapply(df, is.numeric)][,-1]
#define function to calculate summary statistics for each column
calculate_summary <- function(col) {
#calculate mean, median, range, standard deviation, skewness, kurtosis
mean_val <- mean(col, na.rm = TRUE)
median_val <- median(col, na.rm = TRUE)
range_val <- max(col, na.rm = TRUE) - min(col, na.rm = TRUE)
std_dev <- sd(col, na.rm = TRUE)
skew_val <- skew(col, na.rm = TRUE)
kurt_val <- kurtosi(col, na.rm = TRUE)
n_val <- length(col)
min_val <- min(col, na.rm = TRUE)
max_val <- max(col, na.rm = TRUE)
class_val <- class(col)
#return summary stats as a list
return(list(mean = mean_val, median = median_val, range = range_val,
standard_deviation = std_dev, skewness = skew_val, kurtosis = kurt_val,
n = n_val, min = min_val, max = max_val, class = class_val))
}
#apply the function to each numeric column
summary_list <- lapply(numeric_cols, calculate_summary)
#combine results into a dataframe
summary_df <- do.call(rbind, summary_list)
colnames(summary_df) <- c("Mean", "Median", "Range", "Standard_Deviation",
"Skewness", "Kurtosis", "N", "Min", "Max", "Class")
#return the dataframe
return(summary_df)
}
#calculate stats for our dataframe
stats <- summary_stats(df)
#looking at these stats, I think these would be good outcome variables to look at
#psar means postive autonomy and relatedness behaviors (pos AR)
#ngar means negative autonomy and relatedness behaviors (neg AR)
#m means mom and d means dad so these variables represent the pos & neg relatedness behaviors exhibited by participants in interactions with their mothers and fathers
outcome_variables <- c("psard1y", "psarm1y", "ngard1y", "ngarm1y")
#print summary statistics for outcome variables
print(stats[outcome_variables, ])
## Mean Median Range Standard_Deviation Skewness Kurtosis
## psard1y 1.984042 2.041667 2.75 0.5606477 -0.5095955 0.2743636
## psarm1y 1.826142 1.875 3.0625 0.5412087 -0.522489 0.7611735
## ngard1y 0.4093986 0.3 1.95 0.3833093 1.643096 2.999001
## ngarm1y 0.5730465 0.4833333 2.466667 0.4271031 1.370731 2.549321
## N Min Max Class
## psard1y 184 0.375 3.125 "numeric"
## psarm1y 184 0 3.0625 "numeric"
## ngard1y 184 0 1.95 "numeric"
## ngarm1y 184 0 2.466667 "numeric"
#let's make sure these outcome variables look good
for (variable in outcome_variables) {
qqnorm(df[[variable]], main = paste("Q-Q Plot for", variable))
qqline(df[[variable]], col = "red")
}
#and they all look pretty good
#now we can look at predictor variables using a data driven approach
#don't use the id column
numeric_df <- df[, sapply(df, is.numeric)][,-1]
#create an empty list to store all correlation test results
correlation_results <- list()
#loop through each outcome variable
for (outcome_var in outcome_variables) {
#create an empty list to store significant correlations for this outcome variable
significant_correlations <- list()
#loop through each numeric variable in the dataset
for (var_name in colnames(numeric_df)) {
#perform correlation test
cor_result <- cor.test(df[[outcome_var]], df[[var_name]], method = "pearson")
#check if absolute correlation coefficient is above 0.2
if (abs(cor_result$estimate) > 0.2) {
#save significant correlation results
significant_correlations[[var_name]] <- cor_result
}
}
#save significant correlations for this outcome variable in the main list
correlation_results[[outcome_var]] <- significant_correlations
}
#now that we have our correlation results, since there are a lot we should make functions to more easily each through them
#function to print correlation results for a specific variable
print_correlation_results <- function(variable_name) {
if (variable_name %in% names(correlation_results)) {
cat("Correlation results for", variable_name, ":\n")
for (var_name in names(correlation_results[[variable_name]])) {
cor_result <- correlation_results[[variable_name]][[var_name]]
cat("Outcome variable:", var_name, "\n")
cat("Correlation coefficient:", cor_result$estimate, "\n")
cat("p-value:", cor_result$p.value, "\n\n")
}
} else {
cat("No correlation results found for", variable_name, "\n")
}
}
#function to search for a specific variable in the significant correlations of outcome variables
search_predictor_variables <- function(search_variable) {
# Check if the search variable is blank or invalid
if (is.null(search_variable) || trimws(search_variable) == "") {
stop("Search variable is blank or invalid. Please provide a valid search variable.")
}
for (outcome_var in names(correlation_results)) {
significant_correlations <- correlation_results[[outcome_var]]
for (var_name in names(significant_correlations)) {
if (grepl(search_variable, var_name)) {
cat("Found", search_variable, "in significant correlations for", outcome_var, ":\n")
cor_result <- significant_correlations[[var_name]]
cat("Correlation coefficient:", cor_result$estimate, "\n")
cat("p-value:", cor_result$p.value, "\n\n")
}
}
}
}
#i can now search my results to narrow down which predictor variables I might want to use
print_correlation_results("ngard1y")
## Correlation results for ngard1y :
## Outcome variable: cgesc1y
## Correlation coefficient: -0.260147
## p-value: 0.0127608
##
## Outcome variable: deserv1y
## Correlation coefficient: 0.3018021
## p-value: 0.002805159
##
## Outcome variable: blrsd1y
## Correlation coefficient: 0.7227329
## p-value: 6.411785e-17
##
## Outcome variable: valdd1y
## Correlation coefficient: -0.2772184
## p-value: 0.005977964
##
## Outcome variable: enggd1y
## Correlation coefficient: -0.3035564
## p-value: 0.002506222
##
## Outcome variable: intrd1y
## Correlation coefficient: 0.8107307
## p-value: 8.130501e-24
##
## Outcome variable: psrld1y
## Correlation coefficient: -0.323014
## p-value: 0.00125076
##
## Outcome variable: ngatd1y
## Correlation coefficient: 0.9152126
## p-value: 2.822274e-39
##
## Outcome variable: ngrld1y
## Correlation coefficient: 0.8663944
## p-value: 2.099322e-30
##
## Outcome variable: ngard1y
## Correlation coefficient: 1
## p-value: 0
##
## Outcome variable: blrsm1y
## Correlation coefficient: 0.2594436
## p-value: 0.01464868
##
## Outcome variable: enggm1y
## Correlation coefficient: -0.2250349
## p-value: 0.03503538
##
## Outcome variable: intrm1y
## Correlation coefficient: 0.3568769
## p-value: 0.0006423939
##
## Outcome variable: psrlm1y
## Correlation coefficient: -0.2049782
## p-value: 0.05539476
##
## Outcome variable: ngrlm1y
## Correlation coefficient: 0.350251
## p-value: 0.0008219373
##
## Outcome variable: ngarm1y
## Correlation coefficient: 0.2908107
## p-value: 0.005982939
##
## Outcome variable: ycomprd1y
## Correlation coefficient: 0.2203223
## p-value: 0.03191723
##
## Outcome variable: momvrbm1y
## Correlation coefficient: 0.3087638
## p-value: 0.002207609
##
## Outcome variable: alienm1y
## Correlation coefficient: 0.2128937
## p-value: 0.03729345
##
## Outcome variable: pqgrk1y
## Correlation coefficient: 0.257279
## p-value: 0.01138875
##
## Outcome variable: pqfmk1y
## Correlation coefficient: -0.2140578
## p-value: 0.03936739
##
## Outcome variable: engm1y
## Correlation coefficient: -0.2062556
## p-value: 0.05385903
##
## Outcome variable: cism1y
## Correlation coefficient: -0.232827
## p-value: 0.02903754
##
## Outcome variable: satm1y
## Correlation coefficient: -0.3163596
## p-value: 0.002674995
##
## Outcome variable: mcom1y
## Correlation coefficient: -0.3052191
## p-value: 0.003832752
search_predictor_variables("momvrbm1y")
## Found momvrbm1y in significant correlations for ngard1y :
## Correlation coefficient: 0.3087638
## p-value: 0.002207609
#from these, I'm interested in looking at our participant's report of how much verbal abuse they experience from thier mother (momvrbm1y) and their father (dadvrbd1y) and how much the participant's engage with their mothers (enggm1y) and fathers (enggd1y) when their mothers and fathers exhibit autonomy and relatedness behaviors of their own
Visualizing Data
#now that we've identified all the variables we want to use we can start to vizualize them
#first let's look at the distributions for our outcome variables
#this is where we can start to identify abnormalities in our outcomes variables that we might want to fix like outliers
#first melt the data for plotting
melted_data <- reshape2::melt(df[, outcome_variables])
## No id variables; using all as measure variables
## Warning: attributes are not identical across measure variables; they will be
## dropped
#create violin plot with melted data
ggplot(melted_data, aes(x = variable, y = value)) +
geom_violin(trim = FALSE, aes(fill = variable), show.legend = F) +
geom_boxplot(width = 0.1, aes(fill = variable), show.legend = F) +
scale_fill_brewer(palette = "Set2") +
theme_minimal() +
labs(x = "Outcome Variables", y = "Value") +
ggtitle("Violin Plot of Outcome Variables") +
geom_segment(aes(x = 3, y = 3.75, xend = 3.25, yend = 2.5),
arrow = arrow(length = unit(0.5, "cm"))) +
annotate("text", x = 3, y = 4, label = "Neg AR looks a little different than Pos AR",
color = "red", fontface = "bold", size = 3, family = "serif") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
axis.title = element_text(size = 12, family = "serif", face = "bold"),
axis.text.x = element_text(size = 10, family = "serif"),
axis.text.y = element_text(size = 10, family = "serif"))
## Warning: Removed 210 rows containing non-finite values (`stat_ydensity()`).
## Warning: Removed 210 rows containing non-finite values (`stat_boxplot()`).
#from this we can see there might be some outliers and we also get a sense that the neg AR variables are distributed a little differently than the pos AR ones
#let's take a closer look at our outcome variables by making histograms for each one
ngarm1y <- df$ngarm1y[!is.na(df$ngarm1y)]
hist(ngarm1y,
main = "Youths Neg A/R to Mom",
xlab = "Value",
ylab = "Frequency",
col = "orange",
border = "black"
)
ngard1y <- df$ngard1y[!is.na(df$ngard1y)]
hist(ngard1y,
main = "Youths Neg A/R to Dad",
xlab = "Value",
ylab = "Frequency",
col = "orange",
border = "black"
)
psarm1y <- df$psarm1y[!is.na(df$psarm1y)]
hist(psarm1y,
main = "Youths Pos A/R to Mom",
xlab = "Value",
ylab = "Frequency",
col = "orange",
border = "black"
)
psard1y <- df$psard1y[!is.na(df$psard1y)]
hist(psard1y,
main = "Youths Pos A/R to Dad",
xlab = "Value",
ylab = "Frequency",
col = "orange",
border = "black"
)
#now we can explore our some of our correlations look like visualized
ggplot(df, aes(x = momvrbm1y, y = ngarm1y)) +
geom_jitter(color = "#232D4B", alpha = 0.5) + #add points for each data point
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") + #add regression line
labs(x = "Verb Agg from Mom", y = "Child Pos A/R",
title = "Verbal Abuse from Mom in Adolescence vs Neg Autonomy and Relatedness with Mom") +
theme_classic() +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.subtitle = element_text(size = 12, face = "bold", color = "#E57200", hjust = 0.02, family = "serif"),
axis.title = element_text(size = 12, family = "serif", face = "bold"),
axis.text.x = element_text(size = 10, family = "serif"),
axis.text.y = element_text(size = 10, family = "serif"),
panel.grid.minor.y = element_blank(), #get rid of gridlines
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank())
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 19 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 19 rows containing missing values (`geom_point()`).
#we can look at these same scatter plots for all our outcome variables with all our predictor variables
#scatter plots for ngard1y
plot1 <- ggplot(df, aes(x = momvrbm1y, y = ngard1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Mom Verb Abuse", y = "Child Neg A/R to Dad")
plot2 <- ggplot(df, aes(x = dadvrbd1y, y = ngard1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Dad Verb Abuse", y = "Child Neg A/R to Dad")
plot3 <- ggplot(df, aes(x = enggd1y, y = ngard1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Child Eng w/ Dad A/R", y = "Child Neg A/R to Dad")
plot4 <- ggplot(df, aes(x = enggm1y, y = ngard1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Child Eng w/ Mom A/R", y = "Child Neg A/R to Dad")
#scatter plots for ngarm1y
plot5 <- ggplot(df, aes(x = momvrbm1y, y = ngarm1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Mom Verb Abuse", y = "Child Neg A/R to Mom")
plot6 <- ggplot(df, aes(x = dadvrbd1y, y = ngarm1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Dad Verb Abuse", y = "Child Neg A/R to Mom")
plot7 <- ggplot(df, aes(x = enggd1y, y = ngarm1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Child Eng w/ Dad A/R", y = "Child Neg A/R to Mom")
plot8 <- ggplot(df, aes(x = enggm1y, y = ngarm1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Child Eng w/ Mom A/R", y = "Child Neg A/R to Mom")
#arrange neg AR plots together
grid.arrange(plot1, plot2, plot3, plot4, plot5, plot6, plot7, plot8, ncol = 2)
#scatter plots for psard1y
plot1 <- ggplot(df, aes(x = momvrbm1y, y = psard1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Mom Verb Abuse", y = "Child Pos A/R to Dad")
plot2 <- ggplot(df, aes(x = dadvrbd1y, y = psard1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Dad Verb Abuse", y = "Child Pos A/R to Dad")
plot3 <- ggplot(df, aes(x = enggd1y, y = psard1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Child Eng w/ Dad A/R", y = "Child Pos A/R to Dad")
plot4 <- ggplot(df, aes(x = enggm1y, y = psard1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Child Eng w/ Mom A/R", y = "Child Pos A/R to Dad")
# Scatter plots for psarm1y
plot5 <- ggplot(df, aes(x = momvrbm1y, y = psarm1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Mom Verb Abuse", y = "Child Pos A/R to Mom")
plot6 <- ggplot(df, aes(x = dadvrbd1y, y = psarm1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Dad Verb Abuse", y = "Child Pos A/R to Mom")
plot7 <- ggplot(df, aes(x = enggd1y, y = psarm1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Child Eng w/ Dad A/R", y = "Child Pos A/R to Mom")
plot8 <- ggplot(df, aes(x = enggm1y, y = psarm1y)) +
geom_point(color = "#232D4B") +
geom_smooth(method = "lm", se = TRUE, color = "#F84C1E", fill = "#E57200") +
labs(x = "Child Eng w/ Mom A/R", y = "Child Pos A/R to Mom")
#arrange pos AR plots together
grid.arrange(plot1, plot2, plot3, plot4, plot5, plot6, plot7, plot8, ncol = 2)
#let's look at the combined verbal abuse from both parents and the combined engagement with parents from participants
#in order to do that we need to standardized the data, then aggregate the variables
#standardize variables
df_standardized <- df %>%
mutate_at(vars(momvrbm1y, dadvrbd1y, enggm1y, enggd1y), scale)
#take the mean and combine ngard1y and ngarm1y into ngarp1y
df_standardized <- df_standardized %>%
mutate(ngarp1y = rowMeans(dplyr::select(., starts_with("ngar")), na.rm = TRUE))
#take the mean and combine dadvrbd1y and momvrbm1y into parvrbp1y
df_standardized <- df_standardized %>%
mutate(parvrbp1y = rowMeans(dplyr::select(., starts_with("dadvrbd1y"), starts_with("momvrbm1y")), na.rm = TRUE))
#take the mean and combine enggm1y and enggd1y into enggp1y
df_standardized <- df_standardized %>%
mutate(enggp1y = rowMeans(dplyr::select(., starts_with("enggm1y"), starts_with("enggd1y")), na.rm = TRUE))
#let's pick a different set of colors with color brewer
colors <- brewer.pal(3, "Set1")
#scatter plot of parvrbp1y and enggp1y correlated with ngarp1y
ggplot(df_standardized, aes(x = parvrbp1y, y = ngarp1y)) +
geom_point(position = position_jitter(width = 0.1), aes(color = "Parvrbp1y"), alpha = 0.5) +
geom_smooth(method = "lm", aes(color = "Parvrbp1y"), se = TRUE) +
geom_point(aes(x = enggp1y, y = ngarp1y, color = "Enggp1y"), position = position_jitter(width = 0.1), alpha = 0.5) +
geom_smooth(aes(x = enggp1y, y = ngarp1y, color = "Enggp1y"), method = "lm", se = TRUE) +
labs(x = "Mean of Parvrbp1y and Enggp1y (Standardized)",
y = "Mean of Ngarp1y (Standardized)",
title = "Correlation between Parvrbp1y and Enggp1y with Ngarp1y") +
theme_minimal() +
scale_color_manual(values = colors,
labels = c("Parvrbp1y", "Enggp1y")) +
guides(color = guide_legend(title = "Variables"))
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite values (`stat_smooth()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 9 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values (`geom_point()`).
## Warning: Removed 9 rows containing missing values (`geom_point()`).
#the relationship looks a little different when we standardize and aggregate the two parents together
#let's do the exact same thing with the pos AR variables
#take the mean and combine psard1y and psarm1y into psarp1y
df_standardized <- df_standardized %>%
mutate(psarp1y = rowMeans(dplyr::select(., starts_with("psar")), na.rm = TRUE))
#scatter plot of parvrbp1y and enggp1y correlated with psarp1y
ggplot(df_standardized, aes(x = parvrbp1y, y = psarp1y)) +
geom_point(position = position_jitter(width = 0.1), aes(color = "Parvrbp1y"), alpha = 0.5) +
geom_smooth(method = "lm", aes(color = "Parvrbp1y"), se = TRUE) +
geom_point(aes(x = enggp1y, y = psarp1y, color = "Enggp1y"), position = position_jitter(width = 0.1), alpha = 0.5) +
geom_smooth(aes(x = enggp1y, y = psarp1y, color = "Enggp1y"), method = "lm", se = TRUE) +
labs(x = "Mean of Parvrbp1y and Enggp1y (Standardized)",
y = "Mean of Psarp1y (Standardized)",
title = "Correlation between Parvrbp1y and Enggp1y with Psarp1y") +
theme_minimal() +
scale_color_manual(values = colors,
labels = c("Parvrbp1y", "Enggp1y")) +
guides(color = guide_legend(title = "Variables"))
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite values (`stat_smooth()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 9 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values (`geom_point()`).
## Warning: Removed 9 rows containing missing values (`geom_point()`).
#also an interesting relationship: strong postive reltionship between parental abuse and pos AR and a weak relationship between participant's engagement with their parents AR and their pos AR
#really the opposite of what we would expect
Analyzing Data
#let's create a sample linear regression model to explore models may look
#mcontm1y = participant reports their mom continues to argue from the Conflict Tatics Scale (Straus, 1979)
#socacc1y = social acceptance score from Harter's Self-Perception Profile for Children
#ysatc1y = participant's satistication in Supportive Behavior Coding System (Allen et al., 2001)
#mwrm1y = participant's warmth in Supportive Behavior Coding System (Allen et al., 2001)
negARMomLM <- lm(ngarm1y ~ mcontm1y + socacc1y + ysatc1y + mwrm1y, data = df)
summary(negARMomLM)
##
## Call:
## lm(formula = ngarm1y ~ mcontm1y + socacc1y + ysatc1y + mwrm1y,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.82799 -0.22538 -0.08284 0.19929 1.48433
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.37675 0.21291 6.466 1.31e-09 ***
## mcontm1y 0.07503 0.03151 2.381 0.0185 *
## socacc1y -0.02565 0.01105 -2.321 0.0216 *
## ysatc1y -0.08036 0.04876 -1.648 0.1014
## mwrm1y -0.12866 0.05044 -2.551 0.0117 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4037 on 151 degrees of freedom
## (28 observations deleted due to missingness)
## Multiple R-squared: 0.159, Adjusted R-squared: 0.1367
## F-statistic: 7.135 on 4 and 151 DF, p-value: 2.741e-05
plot(negARMomLM)
#everything looks pretty good - not much concern about homogeneity of variance or normality
#no fanning of points in the scale-location plot
#no large curvature in the Q-Q plot
# seems to be pretty much a linear relationship in the Residual vs. Fitted plot
# no overly influential points in the Residuals vs. Leverage plot
#we can use our results here as a baseline for models with our other variables
#there are of course a lot of limitations with this dataset so it is likely the variables will violate some assumptions
#and of course a larger dataset would mean more statistical power
#but let's explore/analyze the data
shapiro.test(df$momvrbm1y)
##
## Shapiro-Wilk normality test
##
## data: df$momvrbm1y
## W = 0.77565, p-value = 2.059e-15
shapiro.test(df$dadvrbd1y)
##
## Shapiro-Wilk normality test
##
## data: df$dadvrbd1y
## W = 0.78534, p-value = 6.014e-14
shapiro.test(df$enggd1y)
##
## Shapiro-Wilk normality test
##
## data: df$enggd1y
## W = 0.96891, p-value = 0.02111
shapiro.test(df$enggm1y)
##
## Shapiro-Wilk normality test
##
## data: df$enggm1y
## W = 0.98035, p-value = 0.0186
shapiro.test(df$ngard1y)
##
## Shapiro-Wilk normality test
##
## data: df$ngard1y
## W = 0.83943, p-value = 7.237e-09
shapiro.test(df$ngarm1y)
##
## Shapiro-Wilk normality test
##
## data: df$ngarm1y
## W = 0.8967, p-value = 2.24e-09
shapiro.test(df$psard1y)
##
## Shapiro-Wilk normality test
##
## data: df$psard1y
## W = 0.97552, p-value = 0.06655
shapiro.test(df$psarm1y)
##
## Shapiro-Wilk normality test
##
## data: df$psarm1y
## W = 0.97912, p-value = 0.01315
#only psard1y doesn't violate our assumption of normality
# leveneTest(df$ngard1y, df$ngarm1y)
# leveneTest(df$psard1y, df$psarm1y)
# leveneTest(df$enggd1y, df$enggm1y)
# leveneTest(df$momvrbm1y, df$dadvrbd1y)
#no violations of homogeneity of variance
#I'm interested in looking at potential disparities between mothers and fathers in the prediction of AR behaviors exhibited by participants
#let's look at whether these might be attributable to differences in the relative levels of observed or reported behaviors for mothers versus fathers
#since assumptions were violated let's do Wilcoxon signed-rank tests
wilcox.test(df$psard1y, df$psarm1y, paired = TRUE)
##
## Wilcoxon signed rank test with continuity correction
##
## data: df$psard1y and df$psarm1y
## V = 1892, p-value = 0.6347
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(df$momvrbm1y, df$dadvrbd1y, paired = TRUE) #it seems there might be a difference in the relative levels of reported verbal abuse from mom and dad, which is important to note for our analysis
##
## Wilcoxon signed rank test with continuity correction
##
## data: df$momvrbm1y and df$dadvrbd1y
## V = 2242.5, p-value = 0.022
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(df$enggd1y, df$enggm1y, paired = TRUE)
##
## Wilcoxon signed rank test with continuity correction
##
## data: df$enggd1y and df$enggm1y
## V = 1421, p-value = 0.2014
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(df$ngard1y, df$ngarm1y, paired = TRUE) #and also a potential difference in the relative levels of observed neg AR towards mom and dad from our participants, which is important to note
##
## Wilcoxon signed rank test with continuity correction
##
## data: df$ngard1y and df$ngarm1y
## V = 1240, p-value = 0.006664
## alternative hypothesis: true location shift is not equal to 0
#we might be able to help our variables look a little better by removing outliers, so let's do that (we can also get a sense of all our variables by making a boxplot for each of them at the same time)
variables <- c("ngard1y", "ngarm1y", "psard1y", "psarm1y", "enggd1y", "enggm1y", "momvrbm1y", "dadvrbd1y")
#create boxplots for each
for (var in variables) {
# Create a boxplot
print(
ggplot(df, aes_string(y = var)) +
geom_boxplot() +
labs(title = paste("Boxplot of", var))
)
#potential outliers
outliers <- boxplot(df[[var]], plot = FALSE)$out
#identified outliers
cat("Outliers for", var, ":", outliers, "\n\n")
}
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 87 rows containing non-finite values (`stat_boxplot()`).
## Outliers for ngard1y : 1.7 1.25 1.366667 1.95 1.4 1.4
## Warning: Removed 18 rows containing non-finite values (`stat_boxplot()`).
## Outliers for ngarm1y : 2.2 1.7 1.54 2.466667 1.7
## Warning: Removed 87 rows containing non-finite values (`stat_boxplot()`).
## Outliers for psard1y : 0.6875 0.875 0.375 0.75 0.5 0.875
## Warning: Removed 18 rows containing non-finite values (`stat_boxplot()`).
## Outliers for psarm1y : 0.375 0 0.375 0.2083333
## Warning: Removed 87 rows containing non-finite values (`stat_boxplot()`).
## Outliers for enggd1y : 0.5 0
## Warning: Removed 18 rows containing non-finite values (`stat_boxplot()`).
## Outliers for enggm1y : 0.5 0 0.5 3.75
## Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
## Outliers for momvrbm1y : 14 19 17 20 13 20 13
## Warning: Removed 26 rows containing non-finite values (`stat_boxplot()`).
## Outliers for dadvrbd1y : 15 15 14 13.2 17 17 13 17 25 19
#function to remove outliers from a variable
remove_outliers <- function(data, variable, outliers) {
data %>%
filter(!variable %in% outliers)
}
#now use the remove_outliers function to remove outliers we identify
for (var in variables) {
#same as before identify the outliers, except this time we're going to remove them instead of printing them
outliers <- boxplot(df[[var]], plot = FALSE)$out
#remove outliers
df2 <- remove_outliers(df, var, outliers)
}
#now that we've removed outliers let's build our models
negARMomLM <- lm(ngarm1y ~ enggd1y + enggm1y + momvrbm1y + dadvrbd1y, data = df)
summary(negARMomLM) #consistent with our visualization, it's the part's engagement with their mom's AR that predicts their neg AR, which is interesting
##
## Call:
## lm(formula = ngarm1y ~ enggd1y + enggm1y + momvrbm1y + dadvrbd1y,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.77124 -0.23992 -0.07156 0.20989 1.65365
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.326181 0.187850 7.060 4.67e-10 ***
## enggd1y -0.043677 0.064944 -0.673 0.50311
## enggm1y -0.298433 0.085359 -3.496 0.00076 ***
## momvrbm1y -0.002744 0.012382 -0.222 0.82518
## dadvrbd1y 0.004917 0.010095 0.487 0.62750
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3797 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## Multiple R-squared: 0.2121, Adjusted R-squared: 0.1741
## F-statistic: 5.586 on 4 and 83 DF, p-value: 0.0004953
#create diagnostic plots
plot(negARMomLM)
#looks pretty ok
negARDadLM <- lm(ngard1y ~ enggd1y + enggm1y + momvrbm1y + dadvrbd1y, data = df)
summary(negARDadLM) #nothing too strong but engagement with dad is an ok predictor
##
## Call:
## lm(formula = ngard1y ~ enggd1y + enggm1y + momvrbm1y + dadvrbd1y,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.69053 -0.24060 -0.06377 0.19162 1.05794
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.835884 0.174039 4.803 6.85e-06 ***
## enggd1y -0.158763 0.060169 -2.639 0.009939 **
## enggm1y -0.033019 0.079083 -0.418 0.677371
## momvrbm1y 0.040076 0.011472 3.493 0.000768 ***
## dadvrbd1y -0.018968 0.009352 -2.028 0.045751 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3518 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## Multiple R-squared: 0.237, Adjusted R-squared: 0.2003
## F-statistic: 6.446 on 4 and 83 DF, p-value: 0.0001442
#create diagnostic plots
plot(negARDadLM) #pretty ok, maybe a non-linear relationship, also observation 94 seems to be a little odd
psarMomLM <- lm(psarm1y ~ enggd1y + enggm1y + momvrbm1y + dadvrbd1y, data = df)
summary(psarMomLM) #mom verb abuse is predictive here but it's not strong
##
## Call:
## lm(formula = psarm1y ~ enggd1y + enggm1y + momvrbm1y + dadvrbd1y,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.76339 -0.22397 -0.02404 0.17896 0.85888
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.385708 0.160334 2.406 0.0184 *
## enggd1y 0.015259 0.055431 0.275 0.7838
## enggm1y 0.647156 0.072855 8.883 1.12e-13 ***
## momvrbm1y 0.023185 0.010569 2.194 0.0310 *
## dadvrbd1y -0.004301 0.008616 -0.499 0.6190
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3241 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## Multiple R-squared: 0.5852, Adjusted R-squared: 0.5652
## F-statistic: 29.27 on 4 and 83 DF, p-value: 3.491e-15
#create diagnostic plots
plot(psarMomLM) #mostly pretty good
psarDadLM <- lm(psard1y ~ enggd1y + enggm1y + momvrbm1y + dadvrbd1y, data = df)
summary(psarDadLM) #engagment with dad is incredibly predictive, which was similar for mom
##
## Call:
## lm(formula = psard1y ~ enggd1y + enggm1y + momvrbm1y + dadvrbd1y,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.81017 -0.16072 0.00033 0.14915 0.73370
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.505503 0.157077 3.218 0.00184 **
## enggd1y 0.635306 0.054305 11.699 < 2e-16 ***
## enggm1y -0.035291 0.071376 -0.494 0.62230
## momvrbm1y 0.013310 0.010354 1.285 0.20220
## dadvrbd1y 0.005794 0.008441 0.686 0.49435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3175 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## Multiple R-squared: 0.6885, Adjusted R-squared: 0.6735
## F-statistic: 45.86 on 4 and 83 DF, p-value: < 2.2e-16
#this makes sense because we would probably expect higher engagement with the parent if they're exhibiting greater pos AR
#create diagnostic plots
plot(psarDadLM) #mostly ok but again it's observation 94 that's a little problematic
#let's explore observation 94
#let's look at our means again
sapply(df[, variables], mean, na.rm = T)
## ngard1y ngarm1y psard1y psarm1y enggd1y enggm1y momvrbm1y dadvrbd1y
## 0.4093986 0.5730465 1.9840421 1.8261421 2.3496564 2.1365892 3.3439560 3.7392405
df[94, variables] #it looks like mom's verbal abuse was very high for this participant
## # A tibble: 1 × 8
## ngard1y ngarm1y psard1y psarm1y enggd1y enggm1y momvrbm1y dadvrbd1y
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.95 0.8 1.94 2.25 0.75 1.75 20 3
negARMomLM <- lm(ngarm1y ~ enggd1y + enggm1y + dadvrbd1y, data = df)
summary(negARMomLM)
##
## Call:
## lm(formula = ngarm1y ~ enggd1y + enggm1y + dadvrbd1y, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.77591 -0.24036 -0.07076 0.21466 1.65773
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.320586 0.185089 7.135 3.17e-10 ***
## enggd1y -0.042861 0.064472 -0.665 0.507997
## enggm1y -0.298239 0.084870 -3.514 0.000714 ***
## dadvrbd1y 0.003817 0.008742 0.437 0.663449
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3776 on 84 degrees of freedom
## (96 observations deleted due to missingness)
## Multiple R-squared: 0.2116, Adjusted R-squared: 0.1835
## F-statistic: 7.516 on 3 and 84 DF, p-value: 0.0001626
plot(negARMomLM)
negARDadLM <- lm(ngard1y ~ enggd1y + enggm1y + dadvrbd1y, data = df)
summary(negARDadLM)
##
## Call:
## lm(formula = ngard1y ~ enggd1y + enggm1y + dadvrbd1y, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.72708 -0.22165 -0.06005 0.12332 1.23189
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.917608 0.183600 4.998 3.12e-06 ***
## enggd1y -0.170681 0.063953 -2.669 0.00913 **
## enggm1y -0.035856 0.084187 -0.426 0.67127
## dadvrbd1y -0.002912 0.008671 -0.336 0.73785
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3745 on 84 degrees of freedom
## (96 observations deleted due to missingness)
## Multiple R-squared: 0.1249, Adjusted R-squared: 0.09359
## F-statistic: 3.995 on 3 and 84 DF, p-value: 0.01036
plot(negARDadLM)
psarMomLM <- lm(psarm1y ~ enggd1y + enggm1y + dadvrbd1y, data = df)
summary(psarMomLM)
##
## Call:
## lm(formula = psarm1y ~ enggd1y + enggm1y + dadvrbd1y, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.76854 -0.22509 0.00288 0.15074 0.78684
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.432987 0.162445 2.665 0.00922 **
## enggd1y 0.008364 0.056584 0.148 0.88284
## enggm1y 0.645515 0.074486 8.666 2.79e-13 ***
## dadvrbd1y 0.004988 0.007672 0.650 0.51736
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3314 on 84 degrees of freedom
## (96 observations deleted due to missingness)
## Multiple R-squared: 0.5612, Adjusted R-squared: 0.5455
## F-statistic: 35.8 on 3 and 84 DF, p-value: 5.29e-15
plot(psarMomLM)
psarDadLM <- lm(psard1y ~ enggd1y + enggm1y + dadvrbd1y, data = df)
summary(psarDadLM)
##
## Call:
## lm(formula = psard1y ~ enggd1y + enggm1y + dadvrbd1y, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.76153 -0.17376 -0.00254 0.15026 0.96137
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.53265 0.15626 3.409 0.001 **
## enggd1y 0.63135 0.05443 11.600 <2e-16 ***
## enggm1y -0.03623 0.07165 -0.506 0.614
## dadvrbd1y 0.01113 0.00738 1.508 0.135
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3187 on 84 degrees of freedom
## (96 observations deleted due to missingness)
## Multiple R-squared: 0.6823, Adjusted R-squared: 0.6709
## F-statistic: 60.12 on 3 and 84 DF, p-value: < 2.2e-16
plot(psarDadLM)
#diagnostic plots look a little better without mom's verbal abuse, which makes a lot sense based on our Wilcoxon signed-rank test for parental verbal abuse
describe(df[outcome_variables])
## vars n mean sd median trimmed mad min max range skew kurtosis
## psard1y 1 97 1.98 0.56 2.04 2.01 0.43 0.38 3.12 2.75 -0.51 0.27
## psarm1y 2 166 1.83 0.54 1.88 1.85 0.46 0.00 3.06 3.06 -0.52 0.76
## ngard1y 3 97 0.41 0.38 0.30 0.35 0.30 0.00 1.95 1.95 1.64 3.00
## ngarm1y 4 166 0.57 0.43 0.48 0.52 0.35 0.00 2.47 2.47 1.37 2.55
## se
## psard1y 0.06
## psarm1y 0.04
## ngard1y 0.04
## ngarm1y 0.03
#based on the skew, kurtosi, and the failure of the shapiro wilk's test we should experiment with non-guassian models
#let's look at the density plots for our variables
plot_title <- "Density Plot of ngarm1y"
density_plot <- density(ngarm1y, na.rm = TRUE)
plot(density_plot, main = plot_title, xlab = "Value", ylab = "Density", col = "blue")
plot_title <- "Density Plot of ngard1y"
density_plot <- density(ngard1y, na.rm = TRUE)
plot(density_plot, main = plot_title, xlab = "Value", ylab = "Density", col = "blue")
plot_title <- "Density Plot of psarm1y"
density_plot <- density(psarm1y, na.rm = TRUE)
plot(density_plot, main = plot_title, xlab = "Value", ylab = "Density", col = "blue")
plot_title <- "Density Plot of psard1y"
density_plot <- density(psard1y, na.rm = TRUE)
plot(density_plot, main = plot_title, xlab = "Value", ylab = "Density", col = "blue")
#Fit GLMs with Poisson distribution
summary(glm(ngard1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df, family = poisson(link = "log")))
##
## Call:
## glm(formula = ngard1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y,
## family = poisson(link = "log"), data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.03126 0.71034 -0.044 0.9649
## momvrbm1y 0.05920 0.03556 1.665 0.0959 .
## dadvrbd1y -0.03005 0.03878 -0.775 0.4383
## enggd1y -0.31180 0.25038 -1.245 0.2130
## enggm1y -0.09276 0.35463 -0.262 0.7937
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 28.531 on 87 degrees of freedom
## Residual deviance: 22.373 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## AIC: Inf
##
## Number of Fisher Scoring iterations: 5
summary(glm(ngarm1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df, family = poisson(link = "log")))
##
## Call:
## glm(formula = ngarm1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y,
## family = poisson(link = "log"), data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.723634 0.595806 1.215 0.2245
## momvrbm1y -0.003778 0.042167 -0.090 0.9286
## dadvrbd1y 0.008988 0.032803 0.274 0.7841
## enggd1y -0.062837 0.220614 -0.285 0.7758
## enggm1y -0.548266 0.301086 -1.821 0.0686 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 25.436 on 87 degrees of freedom
## Residual deviance: 19.607 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## AIC: Inf
##
## Number of Fisher Scoring iterations: 5
summary(glm(psard1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df, family = poisson(link = "log")))
##
## Call:
## glm(formula = psard1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y,
## family = poisson(link = "log"), data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.126648 0.371032 -0.341 0.7328
## momvrbm1y 0.006112 0.023926 0.255 0.7984
## dadvrbd1y 0.002710 0.019693 0.138 0.8906
## enggd1y 0.340412 0.127105 2.678 0.0074 **
## enggm1y -0.021233 0.158138 -0.134 0.8932
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 15.1512 on 87 degrees of freedom
## Residual deviance: 5.4733 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4
summary(glm(psarm1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df, family = poisson(link = "log")))
##
## Call:
## glm(formula = psarm1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y,
## family = poisson(link = "log"), data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.150903 0.368251 -0.410 0.6820
## momvrbm1y 0.011860 0.022890 0.518 0.6044
## dadvrbd1y -0.002313 0.019692 -0.117 0.9065
## enggd1y 0.009522 0.124900 0.076 0.9392
## enggm1y 0.329591 0.160769 2.050 0.0404 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 11.3501 on 87 degrees of freedom
## Residual deviance: 5.0636 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4
#none of these are that convincing
#Fit GLMs with Negative Binomial distribution
summary(glm.nb(ngard1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df))
##
## Call:
## glm.nb(formula = ngard1y ~ momvrbm1y + dadvrbd1y + enggd1y +
## enggm1y, data = df, init.theta = 26652.31709, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.03125 0.71035 -0.044 0.965
## momvrbm1y 0.05920 0.03556 1.665 0.096 .
## dadvrbd1y -0.03005 0.03878 -0.775 0.438
## enggd1y -0.31180 0.25038 -1.245 0.213
## enggm1y -0.09276 0.35464 -0.262 0.794
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(26652.32) family taken to be 1)
##
## Null deviance: 71.906 on 87 degrees of freedom
## Residual deviance: 65.748 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## AIC: 134.07
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 26652
## Std. Err.: 456922
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -122.074
summary(glm.nb(ngarm1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df))
##
## Call:
## glm.nb(formula = ngarm1y ~ momvrbm1y + dadvrbd1y + enggd1y +
## enggm1y, data = df, init.theta = 41923.25433, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.723637 0.595813 1.215 0.2245
## momvrbm1y -0.003777 0.042167 -0.090 0.9286
## dadvrbd1y 0.008988 0.032804 0.274 0.7841
## enggd1y -0.062837 0.220615 -0.285 0.7758
## enggm1y -0.548267 0.301088 -1.821 0.0686 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(41923.25) family taken to be 1)
##
## Null deviance: 68.581 on 87 degrees of freedom
## Residual deviance: 62.751 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## AIC: 151.74
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 41923
## Std. Err.: 781223
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -139.744
summary(glm.nb(psard1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df))
##
## Call:
## glm.nb(formula = psard1y ~ momvrbm1y + dadvrbd1y + enggd1y +
## enggm1y, data = df, init.theta = 471624.1803, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.126648 0.371033 -0.341 0.7328
## momvrbm1y 0.006112 0.023926 0.255 0.7984
## dadvrbd1y 0.002710 0.019693 0.138 0.8906
## enggd1y 0.340413 0.127105 2.678 0.0074 **
## enggm1y -0.021233 0.158138 -0.134 0.8932
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(471624.2) family taken to be 1)
##
## Null deviance: 17.9940 on 87 degrees of freedom
## Residual deviance: 8.3162 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## AIC: 243.69
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 471624
## Std. Err.: 13678678
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -231.689
summary(glm.nb(psarm1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df))
##
## Call:
## glm.nb(formula = psarm1y ~ momvrbm1y + dadvrbd1y + enggd1y +
## enggm1y, data = df, init.theta = 523916.2975, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.150903 0.368252 -0.410 0.6820
## momvrbm1y 0.011860 0.022890 0.518 0.6044
## dadvrbd1y -0.002313 0.019692 -0.117 0.9065
## enggd1y 0.009522 0.124900 0.076 0.9392
## enggm1y 0.329591 0.160770 2.050 0.0404 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(523916.3) family taken to be 1)
##
## Null deviance: 12.0857 on 87 degrees of freedom
## Residual deviance: 5.7992 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## AIC: 242.84
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 523916
## Std. Err.: 16145914
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -230.839
#again not too convincing
# Fit GLMs with Gamma distribution
# summary(glm(ngard1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df, family = Gamma(link = "log")))
# summary(glm(ngarm1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df, family = Gamma(link = "log")))
#non-positive values not allowed
summary(glm(psard1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df, family = Gamma(link = "log")))
##
## Call:
## glm(formula = psard1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y,
## family = Gamma(link = "log"), data = df)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2678348 0.0982264 -2.727 0.0078 **
## momvrbm1y 0.0144697 0.0064747 2.235 0.0281 *
## dadvrbd1y -0.0002525 0.0052784 -0.048 0.9620
## enggd1y 0.3971362 0.0339591 11.695 <2e-16 ***
## enggm1y -0.0255382 0.0446338 -0.572 0.5688
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.03942152)
##
## Null deviance: 9.3809 on 87 degrees of freedom
## Residual deviance: 3.5212 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## AIC: 91.498
##
## Number of Fisher Scoring iterations: 6
summary(glm(psarm1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y, data = df, family = Gamma(link = "log")))
##
## Call:
## glm(formula = psarm1y ~ momvrbm1y + dadvrbd1y + enggd1y + enggm1y,
## family = Gamma(link = "log"), data = df)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.227685 0.093324 -2.440 0.0168 *
## momvrbm1y 0.012233 0.006152 1.989 0.0500 .
## dadvrbd1y -0.001525 0.005015 -0.304 0.7618
## enggd1y 0.013339 0.032264 0.413 0.6804
## enggm1y 0.356566 0.042406 8.408 9.97e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.03558497)
##
## Null deviance: 6.5977 on 87 degrees of freedom
## Residual deviance: 3.1378 on 83 degrees of freedom
## (96 observations deleted due to missingness)
## AIC: 80.816
##
## Number of Fisher Scoring iterations: 5
#pos AR with mom model doesn't look too bad
#let's look at potential mediating effects
predictor_variables <- c("enggd1y", "enggm1y", "momvrbm1y", "dadvrbd1y")
#Model 1: Predict mediator variable (momvrbm1y) from predictor variable (enggm1y)
model1 <- lm(momvrbm1y ~ enggm1y, data = df)
#Model 2: Predict response variable (ngarm1y) from mediator variable (momvrbm1y)
model2 <- lm(ngarm1y ~ momvrbm1y, data = df)
#Model 3: Predict response variable (ngarm1y) from both predictor variable (enggm1y) and mediator variable (momvrbm1y)
model3 <- lm(ngarm1y ~ enggm1y + momvrbm1y, data = df)
#Step 2: Calculate indirect effect
#Indirect effect = (effect of enggm1y on momvrbm1y) * (effect of momvrbm1y on ngarm1y)
indirect_effect <- coef(model1)["enggm1y"] * coef(model2)["momvrbm1y"]
#Step 3: Total effect
#Total effect = effect of enggm1y on ngarm1y
total_effect <- coef(model3)["enggm1y"]
#Step 4: Direct effect
#Direct effect = total effect - indirect effect
direct_effect <- total_effect - indirect_effect
#Step 5: Print results
cat("Indirect Effect:", indirect_effect, "\n")
## Indirect Effect: -0.007049813
cat("Direct Effect:", direct_effect, "\n")
## Direct Effect: -0.2001245
cat("Total Effect:", total_effect, "\n")
## Total Effect: -0.2071743
#let's do that for every combination and print all the results
#loop through each combination
for (outcome_var in outcome_variables) {
for (predictor_var in predictor_variables) {
#skip if predictor variable is the same as response variable
if (predictor_var == outcome_var)
next
#skip if mediator variable is the same as response or predictor variable
for (mediator_var in predictor_variables) {
if (mediator_var == outcome_var || mediator_var == predictor_var)
next
#Model 1: Predict mediator variable from predictor variable
model1 <- lm(as.formula(paste(mediator_var, "~", predictor_var)), data = df)
#Model 2: Predict response variable from mediator variable
model2 <- lm(as.formula(paste(outcome_var, "~", mediator_var)), data = df)
#Model 3: Predict response variable from both predictor and mediator variables
model3 <- lm(as.formula(paste(outcome_var, "~", predictor_var, "+", mediator_var)), data = df)
#calculate indirect effect
indirect_effect <- coef(model1)[predictor_var] * coef(model2)[mediator_var]
#calculate total effect
total_effect <- coef(model3)[predictor_var]
#calculate direct effect
direct_effect <- total_effect - indirect_effect
cat("Response Variable:", outcome_var, "\n")
cat("Predictor Variable:", predictor_var, "\n")
cat("Mediator Variable:", mediator_var, "\n")
cat("Indirect Effect:", indirect_effect, "\n")
cat("Direct Effect:", direct_effect, "\n")
cat("Total Effect:", total_effect, "\n\n")
}
}
}
## Response Variable: psard1y
## Predictor Variable: enggd1y
## Mediator Variable: enggm1y
## Indirect Effect: 0.1667258
## Direct Effect: 0.4547169
## Total Effect: 0.6214427
##
## Response Variable: psard1y
## Predictor Variable: enggd1y
## Mediator Variable: momvrbm1y
## Indirect Effect: -0.0004890695
## Direct Effect: 0.6214219
## Total Effect: 0.6209328
##
## Response Variable: psard1y
## Predictor Variable: enggd1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: 0.006181722
## Direct Effect: 0.612708
## Total Effect: 0.6188898
##
## Response Variable: psard1y
## Predictor Variable: enggm1y
## Mediator Variable: enggd1y
## Indirect Effect: 0.4379254
## Direct Effect: -0.4780232
## Total Effect: -0.04009785
##
## Response Variable: psard1y
## Predictor Variable: enggm1y
## Mediator Variable: momvrbm1y
## Indirect Effect: -0.000574397
## Direct Effect: 0.4111979
## Total Effect: 0.4106235
##
## Response Variable: psard1y
## Predictor Variable: enggm1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: 0.0006612858
## Direct Effect: 0.4071195
## Total Effect: 0.4077808
##
## Response Variable: psard1y
## Predictor Variable: momvrbm1y
## Mediator Variable: enggd1y
## Indirect Effect: -0.0190086
## Direct Effect: 0.03896737
## Total Effect: 0.01995877
##
## Response Variable: psard1y
## Predictor Variable: momvrbm1y
## Mediator Variable: enggm1y
## Indirect Effect: -0.008360062
## Direct Effect: 0.01368964
## Total Effect: 0.005329575
##
## Response Variable: psard1y
## Predictor Variable: momvrbm1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: -0.003327276
## Direct Effect: 0.00726696
## Total Effect: 0.003939684
##
## Response Variable: psard1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: enggd1y
## Indirect Effect: -0.01614085
## Direct Effect: 0.02635609
## Total Effect: 0.01021525
##
## Response Variable: psard1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: enggm1y
## Indirect Effect: -0.0007335172
## Direct Effect: 0.001527169
## Total Effect: 0.0007936519
##
## Response Variable: psard1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: momvrbm1y
## Indirect Effect: 0.0002115228
## Direct Effect: -0.006082592
## Total Effect: -0.00587107
##
## Response Variable: psarm1y
## Predictor Variable: enggd1y
## Mediator Variable: enggm1y
## Indirect Effect: 0.2813804
## Direct Effect: -0.2774567
## Total Effect: 0.003923699
##
## Response Variable: psarm1y
## Predictor Variable: enggd1y
## Mediator Variable: momvrbm1y
## Indirect Effect: 0.004873307
## Direct Effect: 0.2764062
## Total Effect: 0.2812795
##
## Response Variable: psarm1y
## Predictor Variable: enggd1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: -0.007136434
## Direct Effect: 0.2774794
## Total Effect: 0.2703429
##
## Response Variable: psarm1y
## Predictor Variable: enggm1y
## Mediator Variable: enggd1y
## Indirect Effect: 0.1925588
## Direct Effect: 0.4512238
## Total Effect: 0.6437826
##
## Response Variable: psarm1y
## Predictor Variable: enggm1y
## Mediator Variable: momvrbm1y
## Indirect Effect: 0.005723548
## Direct Effect: 0.6846418
## Total Effect: 0.6903654
##
## Response Variable: psarm1y
## Predictor Variable: enggm1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: -0.0007634155
## Direct Effect: 0.6780076
## Total Effect: 0.6772441
##
## Response Variable: psarm1y
## Predictor Variable: momvrbm1y
## Mediator Variable: enggd1y
## Indirect Effect: -0.008358209
## Direct Effect: 0.0267773
## Total Effect: 0.01841909
##
## Response Variable: psarm1y
## Predictor Variable: momvrbm1y
## Mediator Variable: enggm1y
## Indirect Effect: -0.01410914
## Direct Effect: 0.02262909
## Total Effect: 0.008519949
##
## Response Variable: psarm1y
## Predictor Variable: momvrbm1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: 0.003841144
## Direct Effect: -0.01598783
## Total Effect: -0.01214669
##
## Response Variable: psarm1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: enggd1y
## Indirect Effect: -0.00709724
## Direct Effect: 0.00970667
## Total Effect: 0.00260943
##
## Response Variable: psarm1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: enggm1y
## Indirect Effect: -0.001237945
## Direct Effect: 0.009611274
## Total Effect: 0.008373329
##
## Response Variable: psarm1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: momvrbm1y
## Indirect Effect: -0.002107707
## Direct Effect: 0.01373475
## Total Effect: 0.01162704
##
## Response Variable: ngard1y
## Predictor Variable: enggd1y
## Mediator Variable: enggm1y
## Indirect Effect: -0.06381295
## Direct Effect: -0.1042752
## Total Effect: -0.1680881
##
## Response Variable: ngard1y
## Predictor Variable: enggd1y
## Mediator Variable: momvrbm1y
## Indirect Effect: -0.02603282
## Direct Effect: -0.1131557
## Total Effect: -0.1391885
##
## Response Variable: ngard1y
## Predictor Variable: enggd1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: -0.0003653347
## Direct Effect: -0.159791
## Total Effect: -0.1601564
##
## Response Variable: ngard1y
## Predictor Variable: enggm1y
## Mediator Variable: enggd1y
## Indirect Effect: -0.1124383
## Direct Effect: 0.07759395
## Total Effect: -0.03484436
##
## Response Variable: ngard1y
## Predictor Variable: enggm1y
## Mediator Variable: momvrbm1y
## Indirect Effect: -0.03057474
## Direct Effect: -0.1038658
## Total Effect: -0.1344406
##
## Response Variable: ngard1y
## Predictor Variable: enggm1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: -3.908145e-05
## Direct Effect: -0.1558527
## Total Effect: -0.1558918
##
## Response Variable: ngard1y
## Predictor Variable: momvrbm1y
## Mediator Variable: enggd1y
## Indirect Effect: 0.0048805
## Direct Effect: 0.02101214
## Total Effect: 0.02589264
##
## Response Variable: ngard1y
## Predictor Variable: momvrbm1y
## Mediator Variable: enggm1y
## Indirect Effect: 0.003199746
## Direct Effect: 0.02815131
## Total Effect: 0.03135106
##
## Response Variable: ngard1y
## Predictor Variable: momvrbm1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: 0.0001966393
## Direct Effect: 0.03876634
## Total Effect: 0.03896298
##
## Response Variable: ngard1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: enggd1y
## Indirect Effect: 0.004144198
## Direct Effect: -0.008024896
## Total Effect: -0.003880698
##
## Response Variable: ngard1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: enggm1y
## Indirect Effect: 0.0002807478
## Direct Effect: -0.0003992416
## Total Effect: -0.0001184938
##
## Response Variable: ngard1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: momvrbm1y
## Indirect Effect: 0.01125921
## Direct Effect: -0.02645013
## Total Effect: -0.01519093
##
## Response Variable: ngarm1y
## Predictor Variable: enggd1y
## Mediator Variable: enggm1y
## Indirect Effect: -0.08563542
## Direct Effect: 0.03937591
## Total Effect: -0.04625951
##
## Response Variable: ngarm1y
## Predictor Variable: enggd1y
## Mediator Variable: momvrbm1y
## Indirect Effect: -0.006002553
## Direct Effect: -0.1620673
## Total Effect: -0.1680698
##
## Response Variable: ngarm1y
## Predictor Variable: enggd1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: -0.01244055
## Direct Effect: -0.151459
## Total Effect: -0.1638995
##
## Response Variable: ngarm1y
## Predictor Variable: enggm1y
## Mediator Variable: enggd1y
## Indirect Effect: -0.1215692
## Direct Effect: -0.1779959
## Total Effect: -0.2995651
##
## Response Variable: ngarm1y
## Predictor Variable: enggm1y
## Mediator Variable: momvrbm1y
## Indirect Effect: -0.007049813
## Direct Effect: -0.2001245
## Total Effect: -0.2071743
##
## Response Variable: ngarm1y
## Predictor Variable: enggm1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: -0.00133082
## Direct Effect: -0.2214238
## Total Effect: -0.2227546
##
## Response Variable: ngarm1y
## Predictor Variable: momvrbm1y
## Mediator Variable: enggd1y
## Indirect Effect: 0.005276834
## Direct Effect: -0.0040523
## Total Effect: 0.001224534
##
## Response Variable: ngarm1y
## Predictor Variable: momvrbm1y
## Mediator Variable: enggm1y
## Indirect Effect: 0.004293981
## Direct Effect: -0.00157706
## Total Effect: 0.002716921
##
## Response Variable: ngarm1y
## Predictor Variable: momvrbm1y
## Mediator Variable: dadvrbd1y
## Indirect Effect: 0.006696052
## Direct Effect: 0.004124963
## Total Effect: 0.01082102
##
## Response Variable: ngarm1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: enggd1y
## Indirect Effect: 0.004480739
## Direct Effect: 0.0004357268
## Total Effect: 0.004916466
##
## Response Variable: ngarm1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: enggm1y
## Indirect Effect: 0.0003767567
## Direct Effect: 0.01169078
## Total Effect: 0.01206753
##
## Response Variable: ngarm1y
## Predictor Variable: dadvrbd1y
## Mediator Variable: momvrbm1y
## Indirect Effect: 0.002596107
## Direct Effect: 0.005886918
## Total Effect: 0.008483024
#potential mediating effects not very strong except for potentially pos AR with dad and neg AR with dad (which actually seems to outweigh the postive direct effect observed for this model)
#finally we can also look at the models with gender and poverty level (latent variable for socioeconomic status)
#with the amount of missing data we can use Full information maximum likelihood (FIML) to handle missing data
#let's look at this for all our outcome variables
#the gender variable for the age we're looking at is in all caps for some reason when the other ones aren't so let's quickly rename that column
colnames(data)[colnames(data) == "GENDER1Y"] <- "gender1y"
all_variables <- c(outcome_variables, predictor_variables, "gender1y", "perpov1y")
data_subset <- data[, all_variables]
#standardize the data
standardized_data <- scale(data_subset)
standardized_df <- as.data.frame(standardized_data)
colnames(standardized_df) <- all_variables
#perpov is poverty level
#aggregate as we did before with our predictor variables
standardized_df$parvrbp1y <- rowMeans(standardized_df[c("momvrbm1y", "dadvrbd1y")], na.rm = TRUE)
standardized_df$enggp1y <- rowMeans(standardized_df[c("enggd1y", "enggm1y")], na.rm = TRUE)
formula <- as.formula("ngard1y ~ b1*gender1y + b2*perpov1y + b3*momvrbm1y + b4*dadvrbd1y + b5*enggd1y + b6*enggm1y")
lm_model <- lavaan::sem(formula, data = standardized_data, estimator = "ML", missing = "FIML")
## Warning in lavaan::lavaan(model = formula, data = standardized_data, estimator
## = "ML", : lavaan WARNING: model seems to be a formula; please enclose the model
## syntax between quotes
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: 96 cases were deleted due to missing values in
## exogenous variable(s), while fixed.x = TRUE.
print(summary(lm_model))
## lavaan 0.6.17 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 8
##
## Used Total
## Number of observations 88 184
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## ngard1y ~
## gender1y (b1) 0.033 0.098 0.341 0.733
## perpov1y (b2) 0.186 0.143 1.305 0.192
## momvrbm1y (b3) 0.401 0.127 3.147 0.002
## dadvrbd1y (b4) -0.183 0.115 -1.598 0.110
## enggd1y (b5) -0.360 0.119 -3.022 0.003
## enggm1y (b6) -0.017 0.124 -0.137 0.891
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .ngard1y 0.038 0.121 0.316 0.752
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .ngard1y 0.779 0.117 6.633 0.000
formula <- as.formula("ngarm1y ~ b1*gender1y + b2*perpov1y + b3*momvrbm1y + b4*dadvrbd1y + b5*enggd1y + b6*enggm1y")
lm_model <- lavaan::sem(formula, data = standardized_data, estimator = "ML", missing = "FIML")
## Warning in lavaan::lavaan(model = formula, data = standardized_data, estimator = "ML", : lavaan WARNING: model seems to be a formula; please enclose the model syntax between quotes
## Warning in lavaan::lavaan(model = formula, data = standardized_data, estimator = "ML", : lavaan WARNING: 96 cases were deleted due to missing values in
## exogenous variable(s), while fixed.x = TRUE.
print(summary(lm_model))
## lavaan 0.6.17 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 8
##
## Used Total
## Number of observations 88 184
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## ngarm1y ~
## gender1y (b1) 0.020 0.095 0.214 0.830
## perpov1y (b2) 0.152 0.138 1.101 0.271
## momvrbm1y (b3) -0.064 0.124 -0.517 0.605
## dadvrbd1y (b4) 0.091 0.111 0.817 0.414
## enggd1y (b5) -0.118 0.116 -1.017 0.309
## enggm1y (b6) -0.395 0.121 -3.273 0.001
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .ngarm1y -0.016 0.117 -0.139 0.889
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .ngarm1y 0.735 0.111 6.633 0.000
formula <- as.formula("psard1y ~ b1*gender1y + b2*perpov1y + b3*momvrbm1y + b4*dadvrbd1y + b5*enggd1y + b6*enggm1y")
lm_model <- lavaan:: sem(formula, data = standardized_data, estimator = "ML", missing = "FIML")
## Warning in lavaan::lavaan(model = formula, data = standardized_data, estimator = "ML", : lavaan WARNING: model seems to be a formula; please enclose the model syntax between quotes
## Warning in lavaan::lavaan(model = formula, data = standardized_data, estimator = "ML", : lavaan WARNING: 96 cases were deleted due to missing values in
## exogenous variable(s), while fixed.x = TRUE.
print(summary(lm_model))
## lavaan 0.6.17 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 8
##
## Used Total
## Number of observations 88 184
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## psard1y ~
## gender1y (b1) -0.068 0.060 -1.141 0.254
## perpov1y (b2) 0.117 0.087 1.344 0.179
## momvrbm1y (b3) 0.071 0.078 0.914 0.361
## dadvrbd1y (b4) 0.069 0.070 0.986 0.324
## enggd1y (b5) 0.821 0.073 11.274 0.000
## enggm1y (b6) -0.037 0.076 -0.485 0.628
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .psard1y -0.056 0.074 -0.762 0.446
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .psard1y 0.291 0.044 6.633 0.000
formula <- as.formula("psarm1y ~ b1*gender1y + b2*perpov1y + b3*momvrbm1y + b4*dadvrbd1y + b5*enggd1y + b6*enggm1y")
lm_model <- lavaan:: sem(formula, data = standardized_data, estimator = "ML", missing = "FIML")
## Warning in lavaan::lavaan(model = formula, data = standardized_data, estimator = "ML", : lavaan WARNING: model seems to be a formula; please enclose the model syntax between quotes
## Warning in lavaan::lavaan(model = formula, data = standardized_data, estimator = "ML", : lavaan WARNING: 96 cases were deleted due to missing values in
## exogenous variable(s), while fixed.x = TRUE.
print(summary(lm_model)) #socioeconomic status is actually significant here
## lavaan 0.6.17 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 8
##
## Used Total
## Number of observations 88 184
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## psarm1y ~
## gender1y (b1) 0.103 0.062 1.662 0.097
## perpov1y (b2) 0.214 0.090 2.370 0.018
## momvrbm1y (b3) 0.133 0.081 1.650 0.099
## dadvrbd1y (b4) 0.023 0.073 0.312 0.755
## enggd1y (b5) -0.045 0.075 -0.603 0.546
## enggm1y (b6) 0.778 0.079 9.897 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .psarm1y -0.011 0.076 -0.142 0.887
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .psarm1y 0.312 0.047 6.633 0.000
#we see even stronger estimates when using FIML
#these models also suggest gender and socioecomic status may be associated with our outcome variables (it would be worth replicating this study with a larger sample size to look at it more in depth)