Six-week R, programming and statistics masterclass

print(6 + 1)                              # Addition returns 7.
print(2 * (1 + 3))                      # Parentheses evaluated first; returns 8.
print(5 + 3 * 2)                        # Multiplication precedes addition: 11.

x <- sum(1:10)                          # 1:10 is an inclusive integer sequence; sum gives 55.
print(x / 10)                           # Stored result re-used; expected 5.5.
print(ls())                             # List current objects; not the same as printing values.

course_label <- "PALS practice"          # Character/string is enclosed in quotes.
participants <- 66L                      # L makes an integer literal.
sleep_hours <- c(6, 8, 7, 5)            # c() combines values into a numeric vector.
print(sleep_hours[1])                    # R indices start at 1, not 0.
print(sleep_hours[sleep_hours > 6])      # Logical indexing -> c(8, 7).
print(mean(sleep_hours))                 # Expected 6.5 hours, for illustrative data.

crime <- USArrests                        # Built-in public teaching dataset, no CSV needed.
print(dim(crime))                         # Rows and columns: 50 by 4.
print(names(crime))                       # Murder, Assault, UrbanPop, Rape.
print(head(crime))                        # Display the first six rows.

print(max(c(2, 4, 6)))                   # max() needs a vector or individual numbers.
print(crime$Murder)                       # $ extracts the column named Murder.
print(max(crime$Murder))                  # Expected 17.4 in USArrests.
print(max(crime$Assault))                 # Expected 337 in USArrests.
# Interactive help: type ?max or help("max") in the RStudio Console.

rm(crime)                                # Remove only this named object from Environment.
print(exists("crime"))                   # FALSE: object no longer exists.

question <- "Does sleep predict response time?"  # A question is not a causal conclusion.
rt_demo <- c(620, 610, 590, 605)         # Fictional reaction times in milliseconds.
print(data.frame(trial=seq_along(rt_demo), rt_ms=rt_demo))

# Run install.packages("tidyverse") ONCE in Console if not installed.
# library(tidyverse) loads dplyr, tidyr, ggplot2 etc. from your R library.
pinfo <- read.csv("data/participant-info.csv", stringsAsFactors=FALSE)
responses <- read.csv("data/ahi-cesd.csv", stringsAsFactors=FALSE)
print(c(n_participants=nrow(pinfo), n_responses=nrow(responses)))

all_data <- dplyr::inner_join(pinfo, responses, by=c("id","intervention"))
print(c(rows=nrow(all_data), unique_ids=dplyr::n_distinct(all_data$id)))
stopifnot(nrow(all_data) <= nrow(responses))  # Key matching must not explode rows.

selected_vars <- dplyr::select(all_data, id, intervention, ahiTotal, cesdTotal,
                                sex, age, educ, income, occasion)
print(names(selected_vars))               # Exactly nine specified columns.

# MOODLE CORRECTION: worksheet says 1/2 are pre/post; codebook is authoritative:
# occasion 0 = pretest; 1 = posttest; 2 = one-week follow-up.
pre_post <- dplyr::filter(selected_vars, occasion %in% c(0,1))
print(table(pre_post$occasion))           # Inspect sample size/attrition per occasion.

pre_post_asc <- dplyr::arrange(pre_post, ahiTotal)
print(head(dplyr::select(pre_post_asc, id, occasion, ahiTotal)))

overall_summary <- dplyr::summarise(pre_post,
  n=dplyr::n(), mean_ahi=mean(ahiTotal,na.rm=TRUE),
  mean_cesd=mean(cesdTotal,na.rm=TRUE))
print(overall_summary)

occasion_summary <- pre_post_asc |>
  dplyr::group_by(occasion) |>
  dplyr::summarise(n=dplyr::n(), mean_ahi=mean(ahiTotal),
                   mean_cesd=mean(cesdTotal), .groups="drop")
print(occasion_summary)

pre_post <- dplyr::mutate(pre_post, sadness=cesdTotal-ahiTotal)
print(head(dplyr::select(pre_post, sadness, cesdTotal, ahiTotal)))
# DO NOT interpret "sadness" clinically: AHI/CES-D scales have different constructs/ranges.

# MOODLE CORRECTION: worksheet labels conflict with source codebook.
# Correct coding in code-book.pdf: 1=Female, 2=Male.
pre_post <- dplyr::mutate(pre_post,
    sex_label=factor(sex,levels=c(1,2),labels=c("Female","Male")))
print(table(pre_post$sex_label, useNA="ifany"))

paired_pre <- dplyr::filter(all_data, occasion==0) |>
 dplyr::select(id,intervention,ahi_pre=ahiTotal)
paired_post <- dplyr::filter(all_data, occasion==1) |>
 dplyr::select(id,intervention,ahi_post=ahiTotal)
paired_data <- dplyr::inner_join(paired_pre,paired_post,by=c("id","intervention")) |>
 dplyr::mutate(change_ahi=ahi_post-ahi_pre)
print(dplyr::summarise(paired_data, paired_n=dplyr::n(), mean_change=mean(change_ahi)))

qformats <- read.csv("data/Qformats.csv", stringsAsFactors=FALSE)
scoring <- read.csv("data/Scoring.csv", stringsAsFactors=FALSE)
aq_responses <- read.csv("data/AQ.csv", stringsAsFactors=FALSE)
print(c(participants=nrow(aq_responses),questions=nrow(qformats),rules=nrow(scoring)))

aq_long <- tidyr::pivot_longer(aq_responses,
  cols=dplyr::all_of(paste0("Q",1:10)),names_to="Question",values_to="Response")
print(dim(aq_long))                       # Real source: 66*10=660 rows; public demo differs.

aq_with_format <- dplyr::inner_join(aq_long,qformats,by="Question")
aq_scored <- dplyr::inner_join(aq_with_format,scoring,by=c("QFormat","Response"))
stopifnot(nrow(aq_scored)==nrow(aq_long), !anyNA(aq_scored$Score))
print(head(dplyr::select(aq_scored,Id,Question,Response,QFormat,Score)))

aq_scores <- aq_scored |>
  dplyr::group_by(Id) |>
  dplyr::summarise(AQTotal=sum(Score),.groups="drop")
stopifnot(all(aq_scores$AQTotal >= 0 & aq_scores$AQTotal <= 10))
print(head(aq_scores))

aq_pipe <- aq_responses |>
 tidyr::pivot_longer(cols=dplyr::all_of(paste0("Q",1:10)),
                     names_to="Question",values_to="Response") |>
 dplyr::inner_join(qformats,by="Question") |>
 dplyr::inner_join(scoring,by=c("QFormat","Response")) |>
 dplyr::group_by(Id) |>
 dplyr::summarise(AQTotal=sum(Score),.groups="drop")
stopifnot(identical(aq_scores,aq_pipe))   # Native |> parallels tidyverse %>% for these operations.

aq_summary <- dplyr::summarise(aq_scores,
  n=dplyr::n(), mean=mean(AQTotal), sd=sd(AQTotal),
  min=min(AQTotal), max=max(AQTotal), above6=sum(AQTotal>6))
print(aq_summary)
print(dplyr::filter(aq_scores,AQTotal>6))

p_aq <- ggplot2::ggplot(aq_scores,ggplot2::aes(x=AQTotal)) +
 ggplot2::geom_histogram(binwidth=1,boundary=-0.5,colour="white") +
 ggplot2::scale_x_continuous(breaks=0:10) + ggplot2::theme_classic() +
 ggplot2::labs(x="AQ-10 score",y="Number of participants",title="AQ-10 score distribution")
print(p_aq)

variance_aq <- var(aq_scores$AQTotal)    # Sample variance (n-1 denominator).
se_aq <- sd(aq_scores$AQTotal)/sqrt(nrow(aq_scores))
print(c(sample_variance=variance_aq, standard_error=se_aq))

p_hist <- ggplot2::ggplot(all_data,ggplot2::aes(x=ahiTotal)) +
 ggplot2::geom_histogram(binwidth=5,colour="white") +
 ggplot2::labs(x="AHI total",y="Count",title="AHI distribution")
print(p_hist)
print(ggplot2::ggplot(all_data,ggplot2::aes(x=ahiTotal)) +
 ggplot2::geom_histogram(binwidth=1,colour="white") +
 ggplot2::labs(x="AHI total",y="Count",title="AHI distribution: binwidth 1"))

all_data <- dplyr::mutate(all_data,
 sex_label=factor(sex,levels=c(1,2),labels=c("Female","Male")))
p_scatter <- ggplot2::ggplot(all_data,
 ggplot2::aes(x=ahiTotal,y=cesdTotal,colour=sex_label)) +
 ggplot2::geom_point(alpha=0.48) +
 ggplot2::labs(x="AHI total",y="CES-D total",colour="Sex")
print(p_scatter)

p_box <- ggplot2::ggplot(all_data,ggplot2::aes(x=factor(occasion),y=ahiTotal)) +
 ggplot2::geom_boxplot(fill="#7cccd2",outlier.alpha=.4) +
 ggplot2::labs(x="Occasion code",y="AHI total")
print(p_box)
p_violin <- ggplot2::ggplot(all_data,ggplot2::aes(x=factor(occasion),y=ahiTotal)) +
 ggplot2::geom_violin(fill="#c3aaf5",alpha=.55) +
 ggplot2::geom_boxplot(width=.13,outlier.shape=NA) +
 ggplot2::geom_jitter(width=.1,alpha=.15,size=.45)
print(p_violin)

intervention_labels <- c("Signature strengths","Three good things",
                         "Gratitude visit","Early memories (control)")
all_data <- dplyr::mutate(all_data, intervention_label=factor(intervention,
 levels=1:4,labels=intervention_labels))
means_occasion <- all_data |>
 dplyr::group_by(occasion) |>
 dplyr::summarise(mean_ahi=mean(ahiTotal),.groups="drop")
print(ggplot2::ggplot(means_occasion,ggplot2::aes(x=factor(occasion),y=mean_ahi))+
 ggplot2::geom_col(fill="#9fd8c2"))
means_groups <- all_data |>
 dplyr::group_by(occasion,intervention_label) |>
 dplyr::summarise(n=dplyr::n(),mean_ahi=mean(ahiTotal),sd_ahi=sd(ahiTotal),
 .groups="drop")
print(ggplot2::ggplot(means_groups,ggplot2::aes(x=factor(occasion),y=mean_ahi,
 fill=intervention_label)) + ggplot2::geom_col(position="dodge"))

p_line <- ggplot2::ggplot(means_groups,
 ggplot2::aes(x=occasion,y=mean_ahi,colour=intervention_label,group=intervention_label)) +
 ggplot2::geom_line() + ggplot2::geom_point() +
 ggplot2::geom_errorbar(ggplot2::aes(ymin=mean_ahi-sd_ahi,ymax=mean_ahi+sd_ahi),
 width=.13,na.rm=TRUE) + ggplot2::coord_cartesian(ylim=c(24,120)) +
 ggplot2::labs(x="Measurement occasion (0=pre, 1=post)",y="Mean AHI ± 1 SD")
print(p_line)

print(p_scatter + ggplot2::theme_classic() +
 ggplot2::labs(title="Happiness and depressive symptoms"))
print(p_scatter + ggplot2::theme_bw() + ggplot2::facet_wrap(~sex_label))
print(ggplot2::ggplot(all_data,ggplot2::aes(x=ahiTotal,y=cesdTotal)) +
 ggplot2::geom_point(colour="purple",size=2,shape=17,alpha=.3) +
 ggplot2::labs(title="Styled scatterplot",x="AHI",y="CES-D"))

cor_ahi_cesd <- cor(all_data$ahiTotal,all_data$cesdTotal,use="complete.obs")
print(c(correlation=cor_ahi_cesd))
atsq <- read.csv("data/atsq_scores.csv",stringsAsFactors=FALSE)
p_background <- read.csv("data/p_info.csv",stringsAsFactors=FALSE)
atsq_joint <- dplyr::inner_join(atsq,p_background,by="id")
print(c(atsq_rows=nrow(atsq),background_rows=nrow(p_background),joined=nrow(atsq_joint)))

rt_demo <- c(620,610,590,605,615,600,625,595,608,612,598,603)
rt_n <- length(rt_demo)
rt_mean <- mean(rt_demo)
rt_sd <- sd(rt_demo)
rt_se <- rt_sd / sqrt(rt_n)
print(c(n=rt_n,mean_ms=rt_mean,sd_ms=rt_sd,se_ms=rt_se))

set.seed(46)                              # Fix pseudo-random seed for reproducibility.
means_n4 <- replicate(1000,mean(sample(rt_demo,4,replace=TRUE)))
means_n10 <- replicate(1000,mean(sample(rt_demo,10,replace=TRUE)))
print(c(sd_n4=sd(means_n4),sd_n10=sd(means_n10)))
bootstrap_means <- replicate(1000,mean(sample(rt_demo,rt_n,replace=TRUE)))
bootstrap_ci <- quantile(bootstrap_means,c(.025,.975))
print(bootstrap_ci)

par(mfrow=c(1,2))
hist(means_n4,main="Resample size 4",xlab="Mean RT (ms)")
hist(means_n10,main="Resample size 10",xlab="Mean RT (ms)")
par(mfrow=c(1,1))

mu0 <- 600                               # Hypothesised population mean, milliseconds.
t_manual <- (rt_mean-mu0)/(rt_sd/sqrt(rt_n))
df_manual <- rt_n-1
two_sided_p <- 2*pt(-abs(t_manual),df=df_manual)
print(c(t=t_manual,df=df_manual,p=two_sided_p))
model_t <- t.test(rt_demo,mu=mu0,alternative="two.sided",conf.level=.95)
print(model_t)

critical_t <- qt(.975,df=rt_n-1)
ci_manual <- rt_mean + c(-1,1)*critical_t*rt_se
cohens_d <- (rt_mean-mu0)/rt_sd
print(c(ci_lower=ci_manual[1],ci_upper=ci_manual[2],cohens_d=cohens_d))
stopifnot(isTRUE(all.equal(as.numeric(ci_manual),as.numeric(model_t$conf.int),tolerance=1e-8)))

fit <- lm(cesdTotal~ahiTotal,data=all_data)
print(summary(fit)$coefficients)
print(c(r=cor(all_data$ahiTotal,all_data$cesdTotal),
        r_squared=summary(fit)$r.squared))
# Assumptions: independence, approximately linear relationship, conditional errors.