Psi Chi R - May 2025

Start Date: 3 May 2025

Report Date: 31 May 2025

Source: Psi Chi R

In May 2024, Psi Chi R conducted a workshop at APS to help everyone learn R! That workshop was a success, and while we will not be at APS this year, we wanted to re-share that workshop as a resource for you! This month’s contest is a revamped version of what we covered in the workshop. I hope you enjoy! Enjoy the remainder of your semesters – summer is on the horizon!

knitr::opts_chunk$set(echo = T,message = F,warning = F)

#setwd("C:/Users/alanh/Documents/R/Psi_Chi_R")

rm(list=ls())

setwd("~/R/Psi_Chi_R")

library(tidyverse)

#total for bottom row
sum_rows = function(x) {
  x = as.data.frame(x)
  sums = sapply(x,function(col) if (is.numeric(col)) sum(col, na.rm = T) else NA)
  sums = as.data.frame(t(sums))
  names(sums) = names(x)
  rbind(x, sums)
}

## right column for total
sum_cols = function(x) {
  x$Total = rowSums(x[sapply(x, is.numeric)], na.rm = T)
  x
}

#dollar format function
dollars = function(x) {
  paste0("$",format(x,big.mark= ",",scientific=F))
}

desc_stats = function(x){
  c(min = min(x,na.rm=T),
    median = median(x,na.rm=T),
    max = max(x,na.rm=T),
    mean = mean(x,na.rm=T),
    sd = sd(x,na.rm=T))
}

#numeric notations
options(scipen=9999)

Clean and EDA

data=read_csv('data.csv')

names(data) = make.names(colnames(data))

SmartEDA::ExpData(data,type=2) %>% 
  arrange(desc(Per_of_Missing))

Level 1: Preparing the data

#   Remove participants who are missing ‘Age’
#   Remove participants who are missing ‘Gender’
data1= data %>% 
  filter(!is.na(Age) & !is.na(Gender))

Test your skills: Remove participants missing ‘Age’ and ‘Gender’ in one step, instead of two.

#   Create the following variable ‘FACEenmesh’ which measures the extent to which a family is enmeshed, by summing the following items: FACES4, FACES10, FACES16,  FACES22, FACES28, FACES34, FACES40

data2 = data1 %>% 
  mutate(FACEenmesh = FACES4+ FACES10+ FACES16+  FACES22+FACES28+ FACES34+FACES40) %>% 
  select(FACEenmesh,everything())

Level 2: Data inspection + Descriptives

Identify participants with values considered outliers for ‘FACEenmesh’.

data2.1 = data2%>% 
  filter(FACEenmesh > 25)

print(data2.1)

## # A tibble: 12 × 300
##    FACEenmesh ParticipantID   Age Gender Countryborn Stateborn Countrycurrent
##         <dbl>         <dbl> <dbl>  <dbl>       <dbl>     <dbl>          <dbl>
##  1         35             5    19      1         187        25            187
##  2         35            22    18      1         187        25            187
##  3         26            23    18      1         187        25            187
##  4         35            43    19      2         187        44            187
##  5         26            60    18      2         187        25            187
##  6         26           202    18      1         187        51            187
##  7         26           288    18      1         187        25            187
##  8         28           343    19      1         187        44            187
##  9         35           349    18      1         187        11            187
## 10         28           408    18      1         187        25            187
## 11         29           445    19      1         187        25            187
## 12         29           463    19      1         187        25            187
## # ℹ 293 more variables: Statecurrent <dbl>, Yearscurrent <chr>,
## #   Ethnicity <chr>, Race <dbl>, OtherText <chr>, Primarychild <dbl>,
## #   Primarynow <dbl>, Samesex <dbl>, Parentdecease <dbl>,
## #   Parentdeceasetext <chr>, Siblings <dbl>, Primmomnow <chr>,
## #   Primdadnow <chr>, Hoursmom <dbl>, Hoursdad <dbl>, SubSES <dbl>,
## #   FatherEd <dbl>, FatherEdText <chr>, MotherEd <dbl>, MotherEdText <chr>,
## #   Religionnow <dbl>, ReligChristother <chr>, ReligOther <chr>, …

Test your skills: Visualize outlier values for ‘FACEenmesh’.

boxplot(data2$FACEenmesh)
title(main="Looking for FACEenmesh outliers")

Remove participants with values considered outliers on ‘FACEenmesh’; Find the mean of ‘FACEenmesh’; Find the standard deviation of ‘FACEenmesh’; Find the median of ‘FACEenmesh’; Find the range of ‘FACEenmesh’

Test your skills: Find the mean, standard deviation, median, and range of FACEenmesh in one step.

data3 = data2 %>% 
  filter(FACEenmesh <= 25)

#range is from 7 to 25
desc_stats(data3$FACEenmesh) %>% 
  print()

##      min   median      max     mean       sd 
##  7.00000 14.00000 25.00000 14.83878  4.23854

Level 3: Inferential + Other Statistics

Is there a statistically significant correlation between the number of siblings someone has and their enmeshment score?

No, there isn’t.

#No, there isn't a significant correlation.
data4 = data3 %>% 
  filter(!is.na(FACEenmesh) & !is.na(Siblings)) %>% 
  select(Siblings,everything())

cor.test(data4$FACEenmesh,data4$Siblings)

## 
##  Pearson's product-moment correlation
## 
## data:  data4$FACEenmesh and data4$Siblings
## t = -0.35127, df = 457, p-value = 0.7256
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.10779408  0.07521062
## sample estimates:
##         cor 
## -0.01642932

data4 %>% 
  ggplot(aes(y=FACEenmesh,x=Siblings))+
  geom_point()+
  geom_smooth()

Test your skills: Is there a statistically significant difference in enmeshment between people who have siblings and those who do not? *Note, you may need to create a new variable to test this

data5 = data4 %>% 
  mutate(have_siblings = case_when(Siblings > 0 ~ "YES",
                              TRUE ~ 'NO')) %>% 
  select(have_siblings,everything())

shapiro.test(data5$FACEenmesh[data5$have_siblings=='YES'])

## 
##  Shapiro-Wilk normality test
## 
## data:  data5$FACEenmesh[data5$have_siblings == "YES"]
## W = 0.97319, p-value = 0.0000004684

shapiro.test(data5$FACEenmesh[data5$have_siblings=='NO'])

## 
##  Shapiro-Wilk normality test
## 
## data:  data5$FACEenmesh[data5$have_siblings == "NO"]
## W = 0.96121, p-value = 0.2795

hist(data5$FACEenmesh[data5$have_siblings=='YES'])

hist(data5$FACEenmesh[data5$have_siblings=='NO'])

wilcox.test(data5$FACEenmesh~data5$have_siblings)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data5$FACEenmesh by data5$have_siblings
## W = 6810, p-value = 0.7654
## alternative hypothesis: true location shift is not equal to 0

#assumptions for Welch t-test

# data5$dev=with(data5, abs(FACEenmesh - ave(FACEenmesh, have_siblings, FUN = median)))
# 
# summary(aov(data5$dev ~data5$have_siblings))
# 
# t.test(data5$FACEenmesh~data5$have_siblings,var.equal = F)

Level 4: Visualization

Graph the relationship between the number of siblings someone has and their enmeshment score

data5 %>% 
  ggplot(aes(y=FACEenmesh,x=Siblings))+
  theme_bw()+
  geom_point()+
  labs(title='Sibling Rivalry? More Like Sibling Enmeshment',subtitle='A Totally Scientific Look')+
  theme(plot.title = element_text(hjust = .5),
        plot.subtitle = element_text(hjust = .5))+
  geom_smooth()

Test your skills: Visualize the difference in enmeshment between people who have siblings and those who do not

data5 %>% 
  ggplot(aes(x=have_siblings,y=FACEenmesh,fill=have_siblings))+
  geom_col()+
  theme_bw()+
  labs(title='"Only Child Sanctuary" vs. "Sibling Scrum"',subtitle= "The Enmeshment Edition",x='Have Siblings')+
  theme(plot.title = element_text(hjust = .5) , plot.subtitle = element_text(hjust = .5))+
  guides(fill ='none')