We only use non-parametric tests
AEQ scores are essentially in ordinal scale
Chi-squared test of independence (nominal and ordinal (except AEQ score))
Kruskal Wallis test and post hoc (pairwise Wilcox test) (AEQ score and all other numeric data)
The effect sizes of all hypothesis tests except the post-hoc test of Chi-squared test are included.
Chi-squared – Cramer’s V or Cramer’s \(\phi\) (\(\phi_c\))
There is a disagreement between the criteria of Cramer’s V https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6107969/table/tbl2/?report=objectonly
Kruskal Wallis Test by Rank – Eta-squared (\(\eta^2\))
Pairwise Wilcoxon Rank Sum Test – (Cohen) \(r = z/\sqrt{n}\); z statistic; n: sample size number of pairs
\(\phi_c\), \(\eta^2\), and \(r\) of pairwise Wilcoxon are incomparable.
The effect size of Wilcoxon test is not explained well in the original source, but this method of calculation is provided by both R and SPSS:
Rosenthal, R. (1991). Applied Social Research Methods: Meta-analytic procedures for social research Thousand Oaks, CA: SAGE Publications Ltd doi: 10.4135/9781412984997
Spearman Correlation – \(\rho\) (Coefficient)
Interpretation of the Pearson’s and Spearman’s correlation coefficients.
Akoglu H. (2018). User’s guide to correlation coefficients. Turkish journal of emergency medicine, 18(3), 91–93. https://doi.org/10.1016/j.tjem.2018.08.001
You can view significant data by limiting p value and effect size (e.g. p <0.05, effect.size \(\in [0.4,\infty]\)). All of the interpretation reports generated are in APA format.
Example_APA_format_Kruskal
# function that converts group notations (P,S,C) to more readable texts
group.legend = function(data){
new = data %>%
mutate(group = case_when(
group == "control" ~ "Control",
group == "single" ~ "20 minuted 1X/day",
TRUE ~ "10 minutes 2X/day"
),
group = factor(group, levels = c("Control","20 minuted 1X/day","10 minutes 2X/day"))) %>%
return(new)
}
# remove leading 0 from tables, strings
rm_0 <- function(id){
return(ifelse(id < 1, sub("^0+", "", id), id))
}
p.format <- function(p, sig = 0.05){
# APA format of p values:
# 1. 3 decimal places and no leading zeros for values less than 1
# 2. use "p < .001"
str <- ifelse(p <= sig,ifelse(p < 0.001,"p < .001",sprintf("p = %s",rm_0(round(p, digits = 3)))), "insignificant")
return(str)
}
# apply change of group names
i = 1
for(data in datasets){
datasets[[i]] <- group.legend(data)
i = i+1
}
# Modify factor Levels
# Define factor levels of 3; this step is useful for future correlation analysis and data display
levels = as.character(c(1:3))
labels = c("Barely concerned", "Somewhat concerned", "Very Concerned")
# The variable "data" is used for 1.1
# Do not remove it until 1.2 Numeric Data
datasets.raw$pre_study <- datasets.raw$pre_study %>%
rename(improve.score = exam.score) %>%
left_join(select(datasets$other, emory.id, group)) %>%
mutate(
# Separate Race, religious affiliation, family affiliation, and types of meditation practiced using ';'
race = strsplit(race,";"),
religion = strsplit(religion,";"),
religion.family = strsplit(religion.family,";"),
meditation.type = strsplit(meditation.type,";"),
major = strsplit(major,";"),
worry = factor(worry, levels = as.character(c(1:5))),
concern.1 = factor(concern.1, levels=levels),
concern.2 = factor(concern.2, levels=levels),
concern.3 = factor(concern.3, levels=levels),
concern.4 = factor(concern.4, levels=levels),
concern.5 = factor(concern.5, levels=levels),
concern.6 = factor(concern.6, levels=levels),
religious.level = factor(religious.level , levels=as.character(c(1:5))),
meditation.freq = factor(meditation.freq, levels=as.character(c(1:4))),
meditation.times = factor(meditation.times, levels=as.character(c(1:4))),
relax = factor(relax, levels=as.character(c(1:5))),
focus = factor(focus, levels=as.character(c(1:5))),
improve.score = factor(improve.score, levels=as.character(c(1:5))),
)
levels(datasets.raw$pre_study$concern.1) <- c("Barely concerned", "Somewhat concerned", "Very Concerned")
levels(datasets.raw$pre_study$concern.2) <- c("Barely concerned", "Somewhat concerned", "Very Concerned")
levels(datasets.raw$pre_study$concern.3) <- c("Barely concerned", "Somewhat concerned", "Very Concerned")
levels(datasets.raw$pre_study$concern.4) <- c("Barely concerned", "Somewhat concerned", "Very Concerned")
levels(datasets.raw$pre_study$concern.5) <- c("Barely concerned", "Somewhat concerned", "Very Concerned")
levels(datasets.raw$pre_study$concern.6) <- c("Barely concerned", "Somewhat concerned", "Very Concerned")
freq <- c("Never", "Weekly", "Monthly", "Other")
times <- c("Never tried it", "Tried it 1-10 times", "Tried it 10-100 times", "Tried it 100-1000 times")
levels(datasets.raw$pre_study$meditation.freq) <- freq
levels(datasets.raw$pre_study$meditation.times) <- times
likert <- c("Stongly disagree", "Disagree", "Neutral", "Agree", "Strongly agree")
levels(datasets.raw$pre_study$relax) <- likert
levels(datasets.raw$pre_study$focus) <- likert
levels(datasets.raw$pre_study$improve.score) <- likert
demographic <- datasets.raw$pre_study %>%
select(1:35,72)
# remove unnecessary data
preserve <- append(preserve, c("group.legend", "demographic"))
rm(list=ls()[! ls() %in% preserve])
Ordinal variables:
Nested variables
You can interact with the following tables.
demographic.cat.1 <- c("gender", "grade","ethnicity","meditation.freq","meditation.times")
concerns = paste0("concern.",c(1:6))
concerns.label = paste("If you were to get a bad grade, how concerned would you be about:",
c("Disappointing your parents",
"Embarrassing yourself amongst friends",
"Embarrassing yourself amongst professors",
"Not getting the job or into the grad school that you hope to",
"Feeling unintelligent",
"Feeling lazy"))
demographic.cat.1 = c(demographic.cat.1, "difficulty", "worry", concerns, "religious.level",
"relax","focus","improve.score")
demographic.cat.1.label = c(demographic.cat.1, "Difficulty of the selected course", "How worried are you about this exam?",
concerns.label,
"How religious do you consider yourself?",
paste("Meditation helps me", c("relax", "focus","improve my exam score")))
demographic.cat.2 = c("race", "major", "religion", "religion.family", "meditation.type")
baseline.vars <- list(demographic.cat.1, demographic.cat.2)
# function to calculate CI
CI.categorical = function(data, x="group", y){ # colname is a string # y is dependent variable, x is the grouping variable
new <- data %>%
group_by_(x,y) %>% # use the 'standard evaluation' counterpart of `group_by`
# ---*** Warning! some subjects left the question `meditation.freq` in blank ***---
filter(!is.na(.data[[y]])) %>% # To refer to column names that are stored as strings, use the `.data` pronoun
summarise(count = n()) %>%
ungroup() %>%
mutate(sum = sum(count),
proportion = count / sum) %>%
group_by_(x,y) %>% # use the 'standard evaluation' counterpart of `group_by`
mutate(lower = prop.test(count, sum)$conf.int[1], # CI.formula(sum, proportion)[1],
upper = prop.test(count, sum)$conf.int[2], # CI.formula(sum,proportion)[2],
) %>%
mutate(independent = x,dependent = y) %>%
rename_(.dots=setNames(x,"independent.value")) %>%
rename_(.dots=setNames(y,"dependent.value")) %>%
# rename_(.dots=setNames("var1","independent")) %>%
# rename_(.dots=setNames("var2","dependent")) %>%
select(independent,independent.value, dependent, dependent.value, proportion, lower, upper) %>% # remove original count of frequences of categorical value and the sum of all instences
arrange_('dependent.value') %>%
ungroup()
return(new)
}
# function for Chi-squared Post-hoc Testing
post.chi <- function(stdres, method = "bonferroni"){ # Post-hoc test of Chi-squared test
chi_square_values <- stdres^2
p_values <- pchisq(chi_square_values, 1, lower.tail = FALSE) # df = 1
adjusted_p_values <- p_values
for (i in 1:nrow(adjusted_p_values)) {
adjusted_p_values[i, ] <- p.adjust(adjusted_p_values[i,
], method = method, n = ncol(adjusted_p_values) *
nrow(adjusted_p_values))
}
p.df <- as.data.frame(adjusted_p_values) %>% rename(p.adj = Freq)
return(p.df$p.adj)
}
# function for Chi-squared Hypothesis Testing
chi <- function(data, x="group", y, bonf = TRUE){# pass the variable name instead of the variable x, y
tbl = table(data[[x]], data[[y]])
result <- chisq.test(data[[x]], data[[y]])
n = sum(tbl)
q = min(ncol(tbl),nrow(tbl))
effect = sqrt((result$statistic)/(n * (q-1))) # parameter is the degree of freedom, calculates Cramer's V
# Post-hoc test
if(!is.na(result$p.value) & result$p.value <= 0.05){
report <<- sprintf("A chi-square test of independence with continuity correction revealed a significant association between %s and %s, $\\chi^2(%d, N = %d) = %.2f$, $%s$, $\\phi_c = %.2f$.\n", x, y, result$parameter, n, result$statistic, p.format(result$p.value), effect)
p <- post.chi(result$stdres)
post.n <- ncol(result$stdres) * nrow(result$stdres)
stdres <<- as.data.frame(result$stdres) %>% # result data frame of post-hoc test
rename(std_residual = Freq) %>%
rename_(.dots=setNames(1,x)) %>%
rename_(.dots=setNames(2,y)) %>%
mutate(p.adj = p,
sig = ifelse(abs(p.adj) <= 0.05, TRUE, FALSE),
report = case_when(
sig == TRUE ~ sprintf("Post-hoc chi-square test of independence with Bonferroni adjustment revealed a significant association between %s: %s and %s: %s, $\\chi^2(%d, N = %d) = %.2f$, $%s$.\n\n",x, .[[1]],y, .[[2]], 1, post.n, std_residual^2, p.format(p.adj)),
# , $\\phi_c = %.2f$.\n don't know how to calculate the effect size
sig == FALSE ~ NA_character_
))
# concatenate report from post-hoc test
for(i in stdres$report){
if(!is.na(i)) report <<- paste0(report, i)# cat(i)
}
}
# result of Chi-squared test, not including post-hoc test
result_vector = c("Chi-squared", x,y, result$statistic, result$parameter, n, result$p.value, effect)
return(result_vector)
}
# data for first test
CI.demographic.cat = data.frame(matrix(ncol = 7, nrow = 0))
chisq.demographic = data.frame(matrix(ncol = 8, nrow = 0))
colnames(chisq.demographic) = c("test.type","independent", "dependent","statistic", "df", "n", "p.value", "effect.size")
# data for post-hoc test
stdres = NULL
chisq.demographic.post <- data.frame(matrix(ncol = 8, nrow = 0))
colnames(chisq.demographic.post) = c("independent", "independent.value", "dependent", "dependent.value", "std_residual", "p.adj", "report")
report.demographic.cat <- ""
# analyze
for(var in demographic.cat.1){
# CI.categorical
result = CI.categorical(demographic, y=var)
CI.demographic.cat = rbind(CI.demographic.cat, result)
# Hypothesis Testing
result = chi(demographic, y=var)
chisq.demographic[nrow(chisq.demographic) + 1,] = result
# post-hoc test
if(TRUE %in% stdres$sig){
stdres <- stdres %>% filter(sig == TRUE) %>%
mutate(independent = names(stdres)[1], dependent = names(stdres)[2]) %>%
rename_(.dots=setNames(names(stdres)[1],"independent.value")) %>%
rename_(.dots=setNames(names(stdres)[2],"dependent.value")) %>%
select(independent, independent.value, dependent, dependent.value, std_residual, p.adj, report)
chisq.demographic.post = rbind(chisq.demographic.post, stdres)
}
# concatenate report
if(!is.na(result[7]) & result[7] <= 0.05){
report.demographic.cat <- paste(report.demographic.cat, report)
}
gc(verbose = FALSE) # free some RAM
}
# Analyze nested variables
for(var in demographic.cat.2){
# CI
tmp = demographic %>% select(group, var) %>%
unnest(var) # unnest the variable
result = CI.categorical(tmp, y=var)
CI.demographic.cat = rbind(CI.demographic.cat, result)
# Hypothesis testing
result = chi(tmp, y=var)
chisq.demographic[nrow(chisq.demographic) + 1,] = result
# Post-hoc testing
if(TRUE %in% stdres$sig){
stdres <- stdres %>% filter(sig == TRUE) %>%
mutate(independent = names(stdres)[1], dependent = names(stdres)[2]) %>%
rename_(.dots=setNames(names(stdres)[1],"independent.value")) %>%
rename_(.dots=setNames(names(stdres)[2],"dependent.value")) %>%
select(independent, independent.value, dependent, dependent.value, std_residual, p.adj, report)
chisq.groupdiff.post = rbind(chisq.groupdiff.post, stdres)
}
# concatenate report
if(!is.na(result[7]) & result[7] <= 0.05){
report.groupdiff.cat <- paste(report.groupdiff.cat, report)
}
gc(verbose = FALSE) # free some RAM
}
# Display
datatable(
CI.demographic.cat, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
datatable(
chisq.demographic, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
# remove unnecessary data
preserve <- append(preserve,c("CI.categorical","p.format","post.chi", "chi", "chisq.demographic","chisq.demographic.post","report.demographic.cat", "prop.tables.demographic","CI.demographic.cat", "baseline.vars"))
rm(list=ls()[! ls() %in% preserve])
prop.tables.demographic <- NULL
chisq.demographic.sig <- chisq.demographic %>% filter(p.value <=0.05)
for(i in 1:nrow(chisq.demographic.sig)){
row <- chisq.demographic.sig[i, ]
data <- demographic %>%
select(row$independent,row$dependent) %>%
# mutate_if(is.factor, as.integer) %>% # can affect sd()
na.omit() # %>% # apply to both variables
# mutate(.data[[row$independent]] = round(.data[[row$independent]]))
if(row$independent %in% c("race", "major", "religion", "religion.family", "meditation.type")){ # these are all part of nominal variables
data = data %>%
unnest(row$independent)
}
tbl <- round(prop.table(table(data)), digits = 2)
print(tbl)
for(row in 1:nrow(tbl)){
for(col in 1:ncol(tbl)){
tbl[row, col] <- rm_0(tbl[row, col])
}
}
prop.tables.demographic<- append(prop.tables.demographic,list(tbl))
}
## meditation.freq
## group Never Weekly Monthly Other
## Control 0.22 0.00 0.10 0.00
## 20 minuted 1X/day 0.22 0.03 0.10 0.00
## 10 minutes 2X/day 0.19 0.10 0.02 0.00
prop.tables.demographic
## [[1]]
## meditation.freq
## group Never Weekly Monthly Other
## Control .22 .1
## 20 minuted 1X/day .22 .03 .1
## 10 minutes 2X/day .19 .1 .02
Interpretation: - A chi-square test of independence with continuity correction revealed a significant association between group and meditation.freq, \(\chi^2(4, N = 58) = 11.33\), \(p = 0.023\).
Post-hoc chi-square test of independence with Bonferroni adjustment revealed a significant association between group: 10 minutes 2X/day and meditation.freq: Weekly, \(\chi^2(1, N = 9) = 8.38\), \(p = 0.034\).
Most feared class and Reason to meditate are not included currently.
# Define functions
CI.numeric <- function(data, x = "group", y){
new = data %>%
group_by_(x) %>%
filter(!is.na(.data[[y]])) %>%
summarise(
mean = t.test(.data[[y]], conf.level = 0.95)$estimate,
lower = t.test(.data[[y]], conf.level = 0.95)$conf.int[1],
lower = case_when(lower < 0 ~ 0, TRUE ~ lower),
upper = t.test(.data[[y]], conf.level = 0.95)$conf.int[2]) %>%
ungroup() %>%
mutate(variable = y) %>%
select_(x, "variable", "mean", "lower", "upper")
return(new)
}
pairwise.wilcox.test2 <- function(x, g, p.adjust.method = p.adjust.methods, paired = FALSE,
...) # Modified according to the stats package as the original function does not provide test statistic
{
p.adjust.method <- match.arg(p.adjust.method)
DNAME <- paste(deparse1(substitute(x)), "and", deparse1(substitute(g)))
g <- factor(g) #
METHOD <- NULL
# unmodified
compare.levels <- function(i, j) {
xi <- x[as.integer(g) == i]
xj <- x[as.integer(g) == j]
if (is.null(METHOD)) {
wt <- wilcox.test(xi, xj, paired = paired, ...)
METHOD <<- wt$method
wt$p.value
}
else wilcox.test(xi, xj, paired = paired, ...)$p.value
}
PVAL <- pairwise.table(compare.levels, levels(g), p.adjust.method)
# modified
PVAL <- as.data.frame(as.table(PVAL)) %>%
filter(!is.na(Freq)) %>%
rename(p.value = Freq)
# The following 2 functions are used to display test statistic
# modified according to original function
compare.levels.stat <- function(i, j) {
xi <- x[as.integer(g) == i]
xj <- x[as.integer(g) == j]
wilcox.test(xi, xj, paired = paired, ...)$statistic
}
pairwise.table.stat <- function(compare.levels, level.names)
{
ix <- setNames(seq_along(level.names), level.names)
stat <- outer(ix[-1L], ix[-length(ix)], function(ivec, jvec) vapply(seq_along(ivec),
function(k) {
i <- ivec[k]
j <- jvec[k]
if (i > j)
compare.levels.stat(i, j)
else NA_real_
}, 0.05))
# il.tri <- lower.tri(pp, TRUE)
# pp[il.tri] <- p.adjust(pp[il.tri], p.adjust.method)
stat
}
STAT <- pairwise.table.stat(compare.levels, levels(g))
STAT <- as.data.frame(as.table(STAT)) %>%
filter(!is.na(Freq)) %>%
rename(statistic = Freq)
ans <- list(method = METHOD, data.name = DNAME, statistic = STAT, p.value = PVAL,
p.adjust.method = p.adjust.method)
class(ans) <- "pairwise.htest2" # self-modified, avoid repetition of class name
ans
}
kruskal_wilcox <- function(data, x = "group", y, sig = 0.05){ # y is the dependent variable (numeric), x is the grouping variable (categorical)
# Kruskal-Wallis rank sum test
data = data %>% select(x,y) %>% na.omit()
kruskal = kruskal_test(reformulate(x,y), data = data) %>% # y~x
rename(test.type = method,
dependent = .y.,
p.value = p) %>%
mutate(independent = x,
effect.size = kruskal_effsize(data, reformulate(x,y))$effsize) %>% # incorporate effect size
select(test.type, independent, dependent, statistic, df, n, p.value, effect.size)
# Multiple pairwise Wilcox unpaired test between groups
if(kruskal$p.value <= sig){
# discarded because problem arises when encountering many ties
# wilcox = pairwise_wilcox_test(data, reformulate(colname2,colname), p.adjust.method = "BH") %>% # y~x
# rename(y = .y.,
# p.value = p,
# x1 = group1,
# x2 = group2) %>%
# mutate(x = colname2,
# test.type = "pairwise Wilcox",
# effect.size = wilcox_effsize(data, reformulate(colname2,colname))$effsize) %>% # incorporate effect size
# select(test.type,x, x1, x2, y, statistic, n1, n2, p.value, p.adj, effect.size)
wilcox = pairwise.wilcox.test2(data[[y]],data[[x]], p.adjust.method = "BH") # y,x
effect = wilcox_effsize(data, reformulate(x,y))
med <- data %>%
group_by(.data[[x]]) %>%
summarise(mdn = median(.data[[y]], na.rm = T))
wilcox = effect %>%
full_join(wilcox$statistic, by = c("group2"= "Var1","group1" = "Var2")) %>%
full_join(wilcox$p.value, by = c("group2"= "Var1","group1" = "Var2")) %>%
left_join(med, by = c("group1" = x)) %>% rename(mdn1 = mdn) %>%
left_join(med, by = c("group2" = x)) %>% rename(mdn2 = mdn) %>%
rename(dependent = .y.,
independent1 = group1,
independent2 = group2,
effect.size = effsize) %>%
mutate(independent = x,
test.type = "pairwise Wilcox") %>%
select(test.type, independent, independent1, independent2, dependent, statistic, n1, n2, mdn1, mdn2, p.value, effect.size)
return(list(kruskal, wilcox))
}
else{
return(list(kruskal))
}
}
kruskal.report <- function(data){
result <- NULL
# x,y,df, n, statistic, p, effsize
for (i in 1:nrow(data)) {
row = data[i,]
# print(row)
if(row$p.value <= 0.05){
str <- sprintf("A Kruskal-Wallis test showed that the mean rank of %s was significantly different between the levels of %s, $H(%d, N = %d) = %.2f$, $%s$, $\\eta^2 = %.2f $.\n", row$dependent,row$independent, row$df, row$n, row$statistic, p.format(row$p.value), row$effect.size)
# print(str)
result <- append(result, str)
}
}
return(result)
# $\phi_c$, $\eta^2$, and $r$
}
wilcox.report <- function(data){
result <- c()
# x,y,df, n, statistic, p, effsize
for (i in 1:nrow(data)) {
row = data[i,]
# print(row)
if(row$p.value <= 0.05){
str.size <- ifelse(row$n1 == row$n2, sprintf("n1 = n2 = %d", row$n1), sprintf("n1 = %d, n2 = %d", row$n1, row$n2))
str <- sprintf("Post-hoc pairwise Wilcoxon rank sum test with Bonferroni adjustment revealed a significant difference in the mean rank of %s between %s: %s ($mdn1 = %.2f$) and %s: %s ($mdn2 = %.2f$), $W(%s) = %.2f$, $%s$, $\\r = %.2f$.\n", row$dependent,row$independent,row$independent1,row$mdn1, row$independent,row$independent2, row$mdn2, str.size, row$statistic, p.format(row$p.value), row$effect.size)
result <- append(result, str)
}
}
return(result)
}
# Define variable names and labels
demographic.num = c("credits","age",
"GPA.previous", "grade.predict")
demographic.num.labels = c("Credits of the current semester","Age", "Expected grade for the most challenging course", "Estimated GPA of the previous semester")
# Define datasets
CI.demographic.num = data.frame(matrix(ncol = 5, nrow = 0))
kruskal.demographic = data.frame(matrix(ncol = 8, nrow = 0))
wilcox.demographic = data.frame(matrix(ncol = 12, nrow = 0))
for(var in demographic.num){
# CI
result = CI.numeric(demographic, y = var)
CI.demographic.num = rbind(CI.demographic.num, result)
# Hypothesis Testing
result = kruskal_wilcox(demographic %>% select(group, var), y = var)
kruskal.demographic = rbind(kruskal.demographic, result[[1]])
if(length(result) > 1){
wilcox.demographic = rbind(wilcox.demographic, result[[2]])
}
gc(verbose = FALSE) # free some RAM
}
baseline.vars <- append(baseline.vars, list(demographic.num))
# CI_func <- CI.numeric
# remove unnecessary data
preserve <- append(preserve,c("pairwise.wilcox.test2","kruskal_wilcox", "kruskal.report", "wilcox.report", "CI.numeric", "CI.demographic.num", "kruskal.demographic", "wilcox.demographic"))
rm(list=ls()[! ls() %in% preserve])
# Display
datatable(
kruskal.demographic, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
# modified mood.6 and mood.7
productivity67 <- datasets.raw$daily %>%
filter(day == 6) %>% select(emory.id, productivity, day) %>%
bind_rows(select(datasets.raw$daily, emory.id, productivity, day) %>% filter(day == 7) ) %>%
left_join(datasets$other %>% select(emory.id, group), by = "emory.id")
# Create Productivity Dataset separately for day 6 and 7
productivity6 <- productivity67 %>%
filter(day == 6) %>%
select(emory.id,group, productivity) %>%
rename(productivity6 = productivity)
productivity7 <- productivity67 %>%
filter(day == 7) %>%
select(emory.id,group, productivity) %>%
rename(productivity7 = productivity)
productivity6$duplicate <- duplicated(productivity6$emory.id)
productivity6$productivity6[which(productivity6$duplicate == TRUE)-1] <- 3.5 # vary between 3 and 4
productivity6 <- productivity6[-which(productivity6$duplicate == TRUE),]
productivity6$duplicate = NULL
datasets.raw$pre_study <- datasets.raw$pre_study %>%
left_join(select(productivity6,-group),by="emory.id") %>%
left_join(select(productivity7,-group),by="emory.id")
productivity67 <- list(productivity6,productivity7)
# Define Variables Names
productivity = c("productivity6", "productivity7")
productivity.label = c("productivity day 6", "productivity day 7")
baseline.vars <- append(baseline.vars, list(productivity))
# Define data sets
CI.productivity <- data.frame(matrix(ncol = 7, nrow = 0))
chisq.productivity = data.frame(matrix(ncol = 8, nrow = 0))
colnames(chisq.productivity) = c("test.type","independent", "dependent","statistic", "df", "n", "p.value", "effect.size")
# analyze
for(var in productivity){
# CI.categorical
i = which(productivity == var)
result = CI.categorical(productivity67[[i]], y=var)
CI.productivity = rbind(CI.productivity, result)
# Hypothesis Testing
result = chi(productivity67[[i]], x="group", y=var)
chisq.productivity[nrow(chisq.productivity) + 1,] = result
gc(verbose = FALSE) # free some RAM
}
# Display
datatable(
CI.productivity, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
datatable(
chisq.productivity, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
# remove unnecessary data
preserve <- append(preserve,c("chisq.productivity", "CI.productivity", "productivity6", "productivity7"))
rm(list=ls()[! ls() %in% preserve])
# Define variable names and labels
vars = c("exam.score", "study.time.total")
labels = c("Exam score", "Hours studying")
baseline.vars <- append(baseline.vars, list(vars))
# Define data sets
CI.grade_hours = data.frame(matrix(ncol = 5, nrow = 0))
kruskal.grade_hours = data.frame(matrix(ncol = 8, nrow = 0))
wilcox.grade_hours = data.frame(matrix(ncol = 12, nrow = 0))
# analyze
for(var in vars){
# CI
result = CI.numeric(datasets$other, y = var)
CI.grade_hours = rbind(CI.grade_hours, result)
# hypothesis testing
result = kruskal_wilcox(datasets$other, y = var)
kruskal.grade_hours = rbind(kruskal.grade_hours, result[[1]])
if(length(result) > 1){
wilcox.grade_hours = rbind(wilcox.grade_hours, result[[2]])
}
gc(verbose = FALSE) # free some RAM
}
# Display
datatable(
CI.grade_hours, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
datatable(
kruskal.grade_hours, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
datatable(
wilcox.grade_hours, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
# Report
cat(kruskal.report(kruskal.grade_hours))
## A Kruskal-Wallis test showed that the mean rank of exam.score was significantly different between the levels of group, $H(2, N = 44) = 6.39$, $p = .041$, $\eta^2 = 0.11 $.
cat(wilcox.report(wilcox.grade_hours))
## Post-hoc pairwise Wilcoxon rank sum test with Bonferroni adjustment revealed a significant difference in the mean rank of exam.score between group: Control ($mdn1 = 81.00$) and group: 20 minuted 1X/day ($mdn2 = 93.00$), $W(n1 = n2 = 15) = 174.00$, $p = .034$, $\r = 0.47$.
# remove unnecessary data
preserve <- append(preserve,c("CI.grade_hours", "kruskal.grade_hours", "wilcox.grade_hours"))
rm(list=ls()[! ls() %in% preserve])
State Measures are AEQs on day 6, 7, and 14. Trait Measure is the difference between AEQ score on Day 21 and before the study began (day 21 - pre study).
CI.AEQ <- function(data){
new = data %>%
group_by(subscale, group) %>%
summarise(
avg = t.test(score, conf.level = 0.95)$estimate,
L = t.test(score, conf.level = 0.95)$conf.int[1],
U = t.test(score, conf.level = 0.95)$conf.int[2]) %>%
ungroup()
return(new)
}
plot.aeq <- function(data, title = "AEQ", x = "Group", y = "Score"){ # plot a dotplot
subscales = unique(data$subscale)
p = data %>%
# filter(subscale == subscales[i]) %>%
ggplot(mapping = aes(x=group, y=score)) +
geom_dotplot(binaxis='y', stackdir='center', alpha = 0.6) +
facet_grid(rows = vars(subscale)) +
stat_summary(fun.data = "mean_cl_normal",
geom = "errorbar",size = 0.6,
width = .3, color = "black") +
stat_summary(fun = "mean", geom = "point", color = "red", size = 2) +
# geom_errorbar(aes(ymin = , ymax = , fill = group), width = 0.1, size = 0.8) +
labs(# title=title,
y =y, x = x # , caption = "Error bars indicate 95% CI"
)
return(p)
}
hypo.aeq <- function(data){
results = vector(mode = "list", length = 2)
kruskal = data.frame(matrix(ncol = 8, nrow = 0))
wilcox = data.frame(matrix(ncol = 12, nrow = 0))
subscales = unique(data$subscale) # get the subscale names by removing duplicates
for(var in subscales){
new = data %>%
filter(subscale == var) %>%
mutate_(.dots=setNames(list(quote(score)),var)) %>% # standard evaluation https://stackoverflow.com/questions/37259581/r-understanding-standard-evaluation-in-mutate
select_('group', var)
result = kruskal_wilcox(new, y = var)
if(length(result) > 1){
kruskal = rbind(kruskal, result[[1]])
wilcox = rbind(wilcox, result[[2]])
}
else{
kruskal = rbind(kruskal, result[[1]])
}
results = list(kruskal, wilcox)
}
return(results)
}
# CI in tables
measure <- c(paste("Total State AEQ, Day", c(6,7,14)), "Trait AEQ change (Day 21 - Day 0)")
AEQ_6.CI <- datasets$AEQ_6 %>% CI.AEQ() %>%
mutate(AEQ.type = "state", measure = measure[1])
AEQ_7.CI <- datasets$AEQ_7 %>% CI.AEQ() %>%
mutate(AEQ.type = "state", measure = measure[2])
AEQ_14.CI <- datasets$AEQ_14 %>% CI.AEQ() %>%
mutate(AEQ.type = "state", measure = measure[3])
difference.CI <- datasets$difference %>%
group_by(subscale, group) %>%
summarise(
avg = t.test(score, conf.level = 0.95)$estimate,
L = t.test(score, conf.level = 0.95)$conf.int[1],
U = t.test(score, conf.level = 0.95)$conf.int[2]) %>%
ungroup() %>%
mutate(AEQ.type = "trait", measure = measure[4])
CI.AEQ <- rbind(AEQ_6.CI, AEQ_7.CI, AEQ_14.CI, difference.CI)
# CI in graphs
# titles = paste("Total State AEQ Score on Day", c(6, 7, 14))
i = 1
AEQ.graphs = list()
for(data in datasets[c(1,2,3)]){
graph = data %>% plot.aeq(title = measure[i])
AEQ.graphs = append(AEQ.graphs, list(graph))
names(AEQ.graphs)[length(AEQ.graphs)] = c(6, 7, 14)[i]
i = i+1
gc(verbose = FALSE) # free some RAM
}
difference.graph <- datasets$difference %>%
plot.aeq(title = measure[4])
AEQ.graphs = append(AEQ.graphs, list(difference.graph))
names(AEQ.graphs)[length(AEQ.graphs)] = "difference"
AEQ.6714 <- ggarrange(AEQ.graphs[[1]],AEQ.graphs[[2]],AEQ.graphs[[3]],
# labels = c(paste("Day", c(6, 7, 14))),
ncol = 1, nrow = 3, # hjust = c(-0.9, -0.9,-0.9),
vjust = 0,
heights = c(0.8,0.8,0.8))
figure <- ggarrange(AEQ.6714,AEQ.graphs[[4]],
labels = c(paste("Day", c(6, 7, 14)), "21-0"),
ncol = 2, hjust = -1.1,
vjust = 0)
figure.AEQ <- annotate_figure(figure,
top = text_grob("AEQ Scores", color = "black", face = "bold", size = 14),
bottom = text_grob("Error bars indicate 95% CI ", color = "black",
hjust = 1, x = 1, face = "italic", size = 10),
)
figure.AEQ
# hypothesis testing
kruskal.AEQ = data.frame(matrix(ncol = 9, nrow = 0))
wilcox.AEQ = data.frame(matrix(ncol = 13, nrow = 0))
i = 1
for(data in datasets[1:4]){
result = hypo.aeq(data)
result[[1]] = result[[1]] %>%
mutate(measure = measure[i])
if(length(result[[2]]) != 0){
result[[2]] = result[[2]] %>%
mutate(measure = measure[i])
}
kruskal.AEQ = rbind(kruskal.AEQ, result[[1]])
wilcox.AEQ = rbind(wilcox.AEQ, result[[2]])
i = i+1
gc(verbose = FALSE) # free some RAM
}
# Display
datatable(
CI.AEQ, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
datatable(
kruskal.AEQ, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
datatable(
wilcox.AEQ, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
# remove unnecessary data
preserve <- append(preserve,c("CI.AEQ", "kruskal.AEQ", "wilcox.AEQ", "figure.AEQ"))
rm(list=ls()[! ls() %in% preserve])
# Report
cat(kruskal.report(kruskal.AEQ))
## A Kruskal-Wallis test showed that the mean rank of anxiety was significantly different between the levels of group, $H(2, N = 42) = 6.40$, $p = .041$, $\eta^2 = 0.11 $.
## A Kruskal-Wallis test showed that the mean rank of pride was significantly different between the levels of group, $H(2, N = 48) = 7.04$, $p = .03$, $\eta^2 = 0.11 $.
cat(wilcox.report(wilcox.AEQ))
## Post-hoc pairwise Wilcoxon rank sum test with Bonferroni adjustment revealed a significant difference in the mean rank of pride between group: Control ($mdn1 = 12.50$) and group: 20 minuted 1X/day ($mdn2 = 15.00$), $W(n1 = 16, n2 = 17) = 209.00$, $p = .026$, $\r = 0.46$.
A Kruskal-Wallis test showed that the mean rank of anxiety was significantly different between the levels of group, \(H(2, N = 42) = 6.40\), \(p = .041\), $^2 = 0.11 $.
A Kruskal-Wallis test showed that the mean rank of pride was significantly different between the levels of group, \(H(2, N = 48) = 7.04\), \(p = .03\), $^2 = 0.11 $.
95% CI for daily mood
group by duration (1~6 or 7~21) and group: 2 * 3 = 6 CI in total
CI.mood <- datasets$mood %>%
group_by(mood.type, group) %>%
summarise(
avg.mean = t.test(mood.avg, conf.level = 0.95)$estimate,
min.mean = t.test(mood.min, conf.level = 0.95)$estimate,
avg.L = t.test(mood.avg, conf.level = 0.95)$conf.int[1],
avg.U = t.test(mood.avg, conf.level = 0.95)$conf.int[2],
min.L = t.test(mood.min, conf.level = 0.95)$conf.int[1],
min.U = t.test(mood.min, conf.level = 0.95)$conf.int[2]) %>%
ungroup()
avg.graph <- datasets$mood %>%
ggplot(mapping = aes(x=mood.type, y=mood.avg)) +
geom_dotplot(binaxis='y', stackdir='center') +
stat_summary(fun.data = "mean_cl_normal",
geom = "errorbar",size = 1,
width = .3, color = "red") +
stat_summary(fun = "mean", geom = "point", color = "red") +
# geom_errorbar(aes(ymin = , ymax = , fill = group), width = 0.1, size = 0.8) +
facet_wrap(~group) +
labs(# title="4.1 Average Daily Mood across Groups",
x ="Day", y = "Mood (Likert Scale)" #, caption = "Error bars indicate 95% CI"
)
min.graph <- datasets$mood %>%
ggplot(mapping = aes(x=mood.type, y=mood.min)) +
geom_dotplot(binaxis='y', stackdir='center') +
stat_summary(fun.data = "mean_cl_normal",
geom = "errorbar",size = 1,
width = .3, color = "red") +
stat_summary(fun = "mean", geom = "point", color = "red") +
# geom_errorbar(aes(ymin = , ymax = , fill = group), width = 0.1, size = 0.8) +
facet_wrap(~group) +
labs( # title="4.2 Minimum Daily Mood across Groups",
x ="Day", y = "Mood (Likert Scale)" # , caption = "Error bars indicate 95% CI"
)
# combine average mood and min mood score into one graph
figure <- ggarrange(avg.graph,min.graph,
labels = c("Average","Minimum"),
ncol = 1, nrow = 2, hjust = c(-1, -0.9), vjust = 0.2)
figure.mood <- annotate_figure(figure,
top = text_grob("Mood Scores", color = "black", face = "bold", size = 14),
bottom = text_grob("Error bars indicate 95% CI ", color = "black",
hjust = 1, x = 1, face = "italic", size = 10),
)
figure.mood
# use plotly to display the gradual change of mood score
p <- datasets.raw$daily %>%
left_join(select(datasets$other, emory.id, group), by = "emory.id") %>%
filter(day > 0) %>%
ggplot(aes(group, mood, frame = day)) +
geom_jitter(position=position_jitter(0.2), cex=1.2) +
stat_summary(fun.data = "mean_cl_normal",
geom = "errorbar",size = 1,
width = .3, color = "red") +
stat_summary(fun = "mean", geom = "point", color = "red") +
labs(title = "Mood Score", x ="Group", y = "Score (Likert Scale)")
figure.mood.ani <- ggplotly(p) %>%
animation_opts(
1000, easing = "elastic", redraw = FALSE
) %>%
animation_slider(
currentvalue = list(prefix = "Day ", font = list(color="black"))
) %>%
layout(annotations = list(x = 1, y = -0.13, text = "Error bars indicate 95% CI",
showarrow = F, xref='paper', yref='paper',
xanchor='right', yanchor='auto', xshift=0, yshift=0,
font=list(size=12, color="gray"))
)
figure.mood.ani
# remove unnecessary data
preserve <- append(preserve,c("CI.mood", "figure.mood", "figure.mood.ani"))
rm(list=ls()[! ls() %in% preserve])
No significant differences in the daily, mean, and minimum mood scores were found across groups.
Since there are many variables (thousands of combinations), insignificant data are excluded from all graphs and tables in this section.
# continuous vs continuous
# continuous <- c("GPA.previous", "grade.predict", "study.time.total")
# mean and min, treat as continuous
# continuous.2 <- c("mood.avg", "mood.min", "productivity.avg")
# discrete <- c("age", "credits", "exam.score")
# list variable names
concerns = paste0("concern.",c(1:6))
ordinal <- c("difficulty", "worry", concerns, "relax", "focus", "improve.score","meditation.freq","meditation.times", "religious.level", "productivity6", "productivity7", "mood.min.1_6", "mood.min.7_21")
nominal <- c("race", "ethnicity", "gender", "grade","ethnicity","major","religion", "religion.family", "meditation.type","course")
# move AEQ subscales from values in `subscale` to column names
pivot_wider_AEQ <- function(data, name){
AEQ.name = paste0("AEQ.",name, ".")
new = select(data, emory.id, subscale, score) %>%
pivot_wider(names_prefix = AEQ.name,
names_from = "subscale",
values_from = "score",
names_sep = ".")
return(new)
}
# prepare data for spearman correlation
spearman.ori <- datasets.raw$pre_study %>%
select(1:35) %>%
# get AEQ scores
left_join(pivot_wider_AEQ(datasets$AEQ_6,6), by = "emory.id") %>%
left_join(pivot_wider_AEQ(datasets$AEQ_7,7), by = "emory.id") %>%
left_join(pivot_wider_AEQ(datasets$AEQ_14,14), by = "emory.id") %>%
left_join(pivot_wider_AEQ(datasets$difference,"difference"), by = "emory.id") %>%
# get mood avg and min scores
left_join(select(datasets$mood, emory.id, mood.type, mood.avg, mood.min) %>%
pivot_wider(names_from = "mood.type",
values_from = c("mood.avg","mood.min"),
names_sep = "."
),
by = "emory.id") %>%
# get other:
left_join(select(datasets$other,-start.date,-group),by = "emory.id") %>%
select(-(1:6)) %>%
select(-meditation.reason, -exam.date) %>%
select(!one_of(nominal)) %>%
mutate_if(is.factor, as.numeric) # the factor can affect the construction of correlation matrix
# function of correlation and contructring correlation matrix plot
# adapted from
# https://towardsdatascience.com/how-to-create-a-correlation-matrix-with-too-many-variables-309cc0c0a57
corr_simple <- function(data, sig=0.05,type = "spearman"){
#run a correlation and drop the insignificant ones
corr <- rcorr(as.matrix(data), type = type)
#prepare to drop duplicates and correlations of 1
corr$P[lower.tri(corr$P,diag=TRUE)] <- NA
#drop perfect correlations
# corr$P[corr$P == 1] <- NA
corr$P <- as.data.frame(as.table(corr$P)) %>% #turn into a 3-column table, important
na.omit() %>% #remove the NA values from above
arrange(Var1,Var2) %>% # arrange the varaibles in alphabetical order
filter(abs(Freq) <= sig) %>% #select significant values
# exclude relationships between AEQ, mood, and concern
filter(!(substr(Var1,1,4) == "AEQ." & substr(Var2,1,4) == "AEQ.")) %>%
filter(!(substr(Var1,1,5) == "mood." & substr(Var2,1,5) == "mood.")) %>%
filter(!(substr(Var1,1,8) == "concern." & substr(Var2,1,8) == "concern."))
#sort by lowest P value
corr$P <- corr$P[order(abs(corr$P$Freq)),]
# remove values from matrix P and n and convert p and n to dataframe
# corr$P <- corr$r %>%
# select(Var1,Var2)%>%
# inner_join(as.data.frame(as.table(corr$P)), by = c("Var1","Var2"))
# corr$n <- corr$r %>%
# select(Var1,Var2)%>%
# inner_join(as.data.frame(as.table(corr$n)) , by = c("Var1","Var2"))
# concatenate r, P, and n into one dataframe
result <- corr$P %>%
rename(p = Freq) %>%
inner_join(as.data.frame(as.table(corr$r)) %>% rename(r = Freq), by = c("Var1","Var2")) %>%
inner_join(as.data.frame(as.table(corr$n)) %>% rename(n = Freq), by = c("Var1","Var2")) %>% # number of pairwise cases
mutate(df = n-2)
# plot the correlation coefficient, create as a global variable
spearman.plt <<- result %>%
select(Var1,Var2,r) %>%
reshape2::acast(Var1~Var2, value.var="r") %>% # shape correlation coefficient into matrix
ggcorrplot(
outline.col = "white", lab = F)
return(result)
}
spearman <- corr_simple(spearman.ori)
m <- list(
l = 0,
r = 0,
b = 100,
t = 0,
pad = 50
)
# These two lists corresponds to values in columns Var1 and Var2. A for loop is used to discard these relationships.
useless1 <- c("AEQ.14.pride","difficulty","worry","AEQ.7.shame", "AEQ.14.pride","AEQ.14.pride",
"difficulty", "meditation.freq",
# not sure
"grade.predict", "AEQ.7.pride", "AEQ.difference.shame")
useless2 <- c("exam.score", "grade.predict","grade.predict","exam.score","mood.min.7~21","mood.min.1~6",
"worry", "meditation.times",
# not sure
"worry", "exam.score", "exam.score")
# remove meaningless pairs
spearman$useless = FALSE
for(i in 1:length(useless1)){
for(j in 1:nrow(spearman)) {
if(spearman[j,]$Var1 == useless1[i] & spearman[j,]$Var2 == useless2[i]){
spearman[j,]$useless = TRUE
}
}
}
spearman <- spearman[-which(spearman$useless == TRUE),] %>%
filter(abs(r) > 0.35) # filter out weak associations
# reconstruct plot
spearman.plt <- spearman %>%
select(Var1,Var2,r) %>%
reshape2::acast(Var1~Var2, value.var="r") %>% # shape correlation coefficient into matrix
ggcorrplot(
outline.col = "white", lab = F)
spearman.plt
# animated plt
ggplotly(spearman.plt) %>% layout(autosize = F, width = 800, height = 700, margin = m)
Spearman Correlation Coefficient in Dataframe and report in APA format (same as the plot)
datatable(
spearman, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
spearman.report <- function(x,y,r,df,p){
str <- sprintf("Spearman’s rank correlation was computed to assess the relationship between %s and %s. There was a %s correlation between the two variables, $r(%d) = %.2f$, $%s$.\n", x,y,ifelse(r > 0, "positive","negative"),df,r,p.format(p))
}
spm.rep <- spearman %>% mutate(report = spearman.report(Var1,Var2,r,df,p))
# cat(spm.rep$report)
All results shown here are filtered by \(p <= 0.05\) and \(|\rho| > 0.35\) . \(\rho\) is correlation coefficient of spearman correlation.
Attitudes toward the effectiveness of meditation are interrelated. There are significant positive associations among the three types of attitudes:
Possible effect of more experience related to meditation:
which may indicate a possible common reason to meditate for those people who believe that meditation can help them relax or mediatation can help those who meditate more frequently and more times.
indicates that people with more experience in mediation may benefit more by having more positive emotions.
predicted grade and meditation frequency. People who meditate more frequently tend to predict a higher score on exam.
Moderate positive relationship related to accuracy of prediction (exam score)
Moderate positive relationship related to study time and productivity
Relationship between age and many AEQ subscales:
worry and AEQ.difference.hopeless (day 21-0) (-)
Examples:
Students who scored a higher level of worriness tend to have less increment in the level of trait hopelessness.
previous GPA and religious level (-)
previous GPA and AEQ.14.shame (-)
AEQ.difference.anger and minimum mood score from day 7 to 21
The relationships between the variable group and other
numeric variables have been analyzed before
confounding.ori <- datasets.raw$pre_study %>%
select(1:35,73:74) %>%
# get AEQ scores
left_join(pivot_wider_AEQ(datasets$AEQ_6,6), by = "emory.id") %>%
left_join(pivot_wider_AEQ(datasets$AEQ_7,7), by = "emory.id") %>%
left_join(pivot_wider_AEQ(datasets$AEQ_14,14), by = "emory.id") %>%
left_join(pivot_wider_AEQ(datasets$difference,"difference"), by = "emory.id") %>%
# get mood avg and min scores
left_join(select(datasets$mood, emory.id, mood.type, mood.avg, mood.min) %>%
mutate(mood.type = case_when(mood.type == "1~6" ~ "1_6", TRUE ~ "7_21")) %>%
pivot_wider(names_from = "mood.type",
values_from = c("mood.avg","mood.min"),
names_sep = "."
),
by = "emory.id") %>%
# get other:
left_join(select(datasets$other,-start.date,-group),by = "emory.id") %>%
# left_join(select(productivity6,-group, -duplicate), by = "emory.id") %>%
# left_join(select(productivity7,-group), by = "emory.id") %>%
select(-(1:6)) %>%
select(-meditation.reason, -exam.date)
numeric <- names(
confounding.ori %>%
select(!one_of(nominal)) %>%
select(!one_of(ordinal))
)
##############################
nominal <- nominal[nominal != "course"] # only delete this variable temporarily
##############################
col.nested = c("race", "major", "religion", "religion.family", "meditation.type")
CI.confounding.cat = data.frame(matrix(ncol = 6, nrow = 0))
chisq.confounding = data.frame(matrix(ncol = 8, nrow = 0))
colnames(chisq.confounding) = c("test.type","independent", "dependent","statistic", "df", "n", "p.value", "effect.size")
chisq.confounding.post <- data.frame(matrix(ncol = 8, nrow = 0))
colnames(chisq.confounding.post) = c("independent", "independent.value", "dependent", "dependent.value", "std_residual", "p.adj", "report")
stdres = NULL
report.confounding.cat <- ""
for(y in ordinal){ # dependent variable
for(x in nominal){ # grouping variable
data <- confounding.ori %>%
select(x,y) %>%
mutate_if(is.factor, as.integer) %>% # can affect sd()
na.omit() # apply to both variables
if(x %in% col.nested){ # these are all part of nominal variables
data = data %>%
unnest(x)
}
# CI.categorical
result = CI.categorical(data, x=x, y=y)
CI.confounding.cat = rbind(CI.confounding.cat, result)
# Hypothesis Testing
result = chi(data, x=x, y=y)
chisq.confounding[nrow(chisq.confounding) + 1,] = result
if(TRUE %in% stdres$sig){
stdres <- stdres %>% filter(sig == TRUE) %>%
mutate(independent = names(stdres)[1], dependent = names(stdres)[2]) %>%
rename_(.dots=setNames(names(stdres)[1],"independent.value")) %>%
rename_(.dots=setNames(names(stdres)[2],"dependent.value")) %>%
select(independent, independent.value, dependent, dependent.value, std_residual, p.adj, report)
chisq.confounding.post = rbind(chisq.confounding.post, stdres)
}
if(!is.na(result[7]) & result[7] <= 0.05){
report.confounding.cat <- paste(report.confounding.cat, report)
}
gc(verbose = FALSE) # free some RAM
}
gc(verbose = FALSE) # free some RAM
}
datatable(
CI.confounding.cat, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
datatable(
chisq.confounding , filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
Filter out insignificant data same as the previous two tables
chisq.sig.confounding <- chisq.confounding %>% filter(p.value <= 0.05)
discard <- chisq.sig.confounding %>% select(independent, dependent)
CI.sig.confounding.cat <- CI.confounding.cat %>%
left_join(chisq.sig.confounding %>% select(independent, dependent, p.value), by = c("independent", "dependent")) %>%
filter(!is.na(p.value))
prop.tables <- NULL
for(i in 1:nrow(chisq.sig.confounding)){
row <- chisq.sig.confounding[i, ]
data <- confounding.ori %>%
select(row$independent,row$dependent) %>%
# mutate_if(is.factor, as.integer) %>% # can affect sd()
na.omit() # %>% # apply to both variables
# mutate(.data[[row$independent]] = round(.data[[row$independent]]))
if(row$independent %in% col.nested){ # these are all part of nominal variables
data = data %>%
unnest(row$independent)
}
tbl <- round(prop.table(table(data)), digits = 2)
for(row in 1:nrow(tbl)){
for(col in 1:ncol(tbl)){
tbl[row, col] <- rm_0(tbl[row, col])
}
}
# print(tbl)
prop.tables<- append(prop.tables,list(tbl))
}
datatable(
CI.sig.confounding.cat, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
datatable(
chisq.sig.confounding, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
datatable(
chisq.confounding.post %>% select(-report), filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
print(prop.tables)
## [[1]]
## concern.2
## grade Barely concerned Somewhat concerned Very Concerned
## Freshman .07 .38 .07
## Sophomore .2 .18 .1
##
## [[2]]
## concern.4
## gender Barely concerned Somewhat concerned Very Concerned
## Female .05 .15 .44
## Male .12 .07 .17
##
## [[3]]
## relax
## meditation.type Stongly disagree Disagree Neutral Agree Strongly agree
## Focused attention .02 .17 .04
## Loving-kindness .02 .03 .02
## Mindfulness .01 .06 .27 .09
## Never tried meditation .07 .01
## Not sure .02 .08 .08
##
## [[4]]
## meditation.freq
## major Never Weekly
## Biology .01
## Anthropology .01
## Biology .06
## Business .01 .01
## Chemistry .04
## Computer Science .01
## Economics .01 .01
## Finance .03 .01
## Human Biology .01
## Human Health .03
## International Studies .01
## Medial Anthropology .01
## Music
## Neuroscience .01
## Neuroscience and Behavioral Biology .17 .01
## Nursing .04
## Physical Therapy .01
## Political Science .01
## Pre-6(Biology Major (Area 6: Pre-Med, Pre-Dental, Pre-Vet)) .01
## Psychology .13 .01
## Quantitative Science
## Spanish
## WGS
## meditation.freq
## major Monthly Other
## Biology
## Anthropology
## Biology .01
## Business .01
## Chemistry
## Computer Science
## Economics .03
## Finance
## Human Biology
## Human Health
## International Studies
## Medial Anthropology
## Music .01
## Neuroscience .01
## Neuroscience and Behavioral Biology .01
## Nursing
## Physical Therapy
## Political Science .03
## Pre-6(Biology Major (Area 6: Pre-Med, Pre-Dental, Pre-Vet))
## Psychology .07
## Quantitative Science .04
## Spanish .01
## WGS .01
##
## [[5]]
## meditation.freq
## meditation.type Never Weekly Monthly Other
## Focused attention .13 .05 .06
## Loving-kindness .02 .02 .04
## Mindfulness .19 .09 .14
## Never tried meditation .08
## Not sure .16 .01
##
## [[6]]
## meditation.times
## gender Never tried it Tried it 1-10 times Tried it 10-100 times
## Female .11 .39 .1
## Male .02 .13 .18
## meditation.times
## gender Tried it 100-1000 times
## Female .05
## Male .02
##
## [[7]]
## meditation.times
## major Never tried it
## Biology
## Anthropology
## Biology .01
## Business
## Chemistry .03
## Computer Science
## Economics
## Finance .01
## Human Biology
## Human Health .01
## International Studies
## Medial Anthropology
## Music
## Neuroscience
## Neuroscience and Behavioral Biology .04
## Nursing
## Physical Therapy
## Political Science
## Pre-6(Biology Major (Area 6: Pre-Med, Pre-Dental, Pre-Vet))
## Psychology .01
## Quantitative Science
## Spanish
## WGS
## meditation.times
## major Tried it 1-10 times
## Biology .01
## Anthropology .01
## Biology .05
## Business .01
## Chemistry .01
## Computer Science
## Economics .01
## Finance .01
## Human Biology .01
## Human Health .01
## International Studies .01
## Medial Anthropology
## Music
## Neuroscience .01
## Neuroscience and Behavioral Biology .14
## Nursing .04
## Physical Therapy
## Political Science .01
## Pre-6(Biology Major (Area 6: Pre-Med, Pre-Dental, Pre-Vet))
## Psychology .11
## Quantitative Science .01
## Spanish
## WGS
## meditation.times
## major Tried it 10-100 times
## Biology
## Anthropology
## Biology
## Business .01
## Chemistry
## Computer Science .01
## Economics .04
## Finance .01
## Human Biology
## Human Health
## International Studies
## Medial Anthropology
## Music .01
## Neuroscience .01
## Neuroscience and Behavioral Biology .03
## Nursing
## Physical Therapy .01
## Political Science .03
## Pre-6(Biology Major (Area 6: Pre-Med, Pre-Dental, Pre-Vet))
## Psychology .08
## Quantitative Science .03
## Spanish
## WGS .01
## meditation.times
## major Tried it 100-1000 times
## Biology
## Anthropology
## Biology
## Business .01
## Chemistry
## Computer Science
## Economics
## Finance
## Human Biology
## Human Health
## International Studies
## Medial Anthropology .01
## Music
## Neuroscience
## Neuroscience and Behavioral Biology
## Nursing
## Physical Therapy
## Political Science
## Pre-6(Biology Major (Area 6: Pre-Med, Pre-Dental, Pre-Vet)) .01
## Psychology
## Quantitative Science .01
## Spanish .01
## WGS
##
## [[8]]
## productivity6
## race 1 2 3 3.5 4 5
## Asian .02 .22 .04 .24 .04
## Black or African American .02 .02 .02
## Other .02 .04
## White .1 .02 .2
##
## [[9]]
## productivity7
## gender 1 2 3 4 5
## Female .02 .04 .23 .33 .02
## Male .12 .02 .1 .08 .02
##
## [[10]]
## mood.min.1_6
## religion.family 1 2 3
## Buddhism .08 .02
## Christianity .04 .19 .11
## Hinduism .02 .06
## Islam .02 .02
## Judaism .04 .02
## None .06 .23 .06
## Other .04
## Prefer not to specify .02
##
## [[11]]
## mood.min.7_21
## race 1 2 3 4
## Asian .1 .33 .13
## Black or African American .02 .04
## Other .02 .02 .02
## White .13 .15 .04
kruskal.confounding = data.frame(matrix(ncol = 8, nrow = 0))
colnames(kruskal.confounding) = c("test.type","independent","dependent","statistic","df","n","p.value","effect.size")
wilcox.confounding = data.frame(matrix(ncol = 12, nrow = 0))
colnames(wilcox.confounding) = c("test.type","independent","independent1","independent2","dependent","statistic","n1","n2","mdn1", "mdn2","p.value","effect.size")
CI.confounding.num = data.frame(matrix(ncol = 6, nrow = 0))
# define new CI
CI.confounding <- function(data, x = "group", y){
new = data %>%
mutate_if(is.factor, as.integer) %>% # can affect sd()
group_by_(x) %>%
filter(!is.na(.data[[y]])) %>%
summarise(
mean = mean(.data[[y]]),
sd = sd(.data[[y]]),
n = n(),
error = qnorm(0.975) - sd/sqrt(n),
lower = mean - error,
lower = case_when(lower < 0 ~ 0, TRUE ~ lower),
upper = mean + error) %>%
ungroup() %>%
mutate(independent = x, dependent = y) %>%
rename(independent.value = x) %>% # need to include a column of grouping variable name since we are dealing with multiple grouping variables instead of one
select(independent,independent.value,dependent, "mean", "lower", "upper")
return(new)
}
# beware of your RAM
for(y in numeric){ # dependent variable
for(x in nominal){ # grouping variable
# subset data set
data = confounding.ori %>%
select(x,y) %>%
mutate_if(is.factor, as.integer) %>% # can affect sd()
na.omit() # apply to both variables
# unnest variables
if(x %in% col.nested){
data = data %>%
unnest(x)
}
# CI
tmp = CI.confounding(data, x, y)
CI.confounding.num = rbind(CI.confounding.num, tmp)
# Hypothesis Testing
tmp = kruskal_wilcox(data, x, y)
kruskal.confounding = rbind(kruskal.confounding, tmp[[1]])
if(length(tmp) > 1){
wilcox.confounding = rbind(wilcox.confounding, tmp[[2]])
}
gc(verbose = FALSE) # free some RAM
}
gc(verbose = FALSE) # free some RAM
}
datatable(
CI.confounding.num, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
datatable(
kruskal.confounding %>% filter(p.value <= 0.05), filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
datatable(
wilcox.confounding %>% filter(p.value <= 0.05), filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
# remove unnecessary data
preserve <- append(preserve,c("pivot_wider_AEQ","spearman", "spearman.plt", "CI.confounding.num", "kruskal.confounding", "wilcox.confounding", "CI.confounding.cat", "chisq.confounding"))
rm(list=ls()[! ls() %in% preserve])
combined_data <- datasets.raw$pre_study %>% # without daily surveys
select(1:35) %>%
# get AEQ scores
left_join(pivot_wider_AEQ(datasets$AEQ_6,6), by = "emory.id") %>%
left_join(pivot_wider_AEQ(datasets$AEQ_7,7), by = "emory.id") %>%
left_join(pivot_wider_AEQ(datasets$AEQ_14,14), by = "emory.id") %>%
left_join(pivot_wider_AEQ(datasets$difference,"difference"), by = "emory.id") %>%
# get mood avg and min scores
left_join(select(datasets$mood, emory.id, mood.type, mood.avg, mood.min) %>%
pivot_wider(names_from = "mood.type",
values_from = c("mood.avg","mood.min"),
names_sep = "."
),
by = "emory.id") %>%
# get other:
left_join(select(datasets$other,emory.id, start.date,group, exam.score, study.time.total),by = "emory.id") %>%
select(-(1:3),-(5:6)) %>%
left_join(select(productivity6,-group), by = "emory.id") %>%
left_join(select(productivity7,-group), by = "emory.id")
# mutate(group = factor(group, levels = c("Control","20 minuted 1X/day","10 minutes 2X/day")))
# select(meditation.reason, exam.date) %>%
# select(!one_of(nominal))
datatable(
combined_data, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
Non-parametric methods of regression are non-linear regression models. Here, generalized additive model is used.
The model generally follows the format of y ~ x1 + x2.
This is a simple analysis of two independent variables and one dependent
variable, namely, AEQ.scores ~ group + baseline.data.
Generalized (Approximate) Cross Validation (GCV or GACV) method is
used.
bs (define smooth terms):
use bs = "fs" for all categorical variables (coerce
them into factors, including nominal variables)
use default setting bs = "tp" for exam score and
hours studying total
From Rdocumentation: “These are low rank isotropic smoothers of any number of covariates. By isotropic is meant that rotation of the covariate co-ordinate system will not change the result of smoothing. By low rank is meant that they have far fewer coefficients than there are data to smooth. They are reduced rank versions of the thin plate splines and use the thin plate spline penalty. They are the default smooth for s terms because there is a defined sense in which they are the optimal smoother of any given basis dimension/rank (Wood, 2003).”
Further model checks should be runned.
gam.check(model)plot(b,pages=1,residuals=TRUE) # show partial residualsplot(b,pages=1,seWithMean=TRUE) # 'with intercept' CIsMore detailed information on this website.
Only significant regression models are included in the output. The data set includes all models.
# AEQs <- c("AEQ_6","AEQ_7","AEQ_14","difference")
AEQs <- names(combined_data)[31:50]
var.types <- c("NA","demographic.cat.1","demographic.cat.2","demographic.num", "productivity", "exam.score_study.time")
regression_data <- combined_data
baseline.vars = append("NA", baseline.vars)
names(baseline.vars) <- var.types
# models <- NULL # 5 dimension
# baseline.vars = baseline.vars[-c(1:11)]
models.df <- data.frame(matrix(ncol = 6, nrow = 0))
colnames(models.df) = c("baseline.var.type","baseline.var","AEQ.subscale","coefficient","smooth.sig")
i = 1
report <- c()
for(baseline.var.type in baseline.vars){
for(baseline.var in baseline.var.type){
# print(baseline.var)
if(var.types[i] == "demographic.cat.2"){
regression_unnest <- regression_data %>%
unnest(baseline.var) %>% na.omit(cols=baseline.var)
}
else regression_unnest <- regression_data %>% na.omit(cols=baseline.var)
if(var.types[i] != "NA" & var.types[i] != "exam.score_study.time"){
regression_unnest[baseline.var] = lapply(regression_unnest[baseline.var], factor)
}
# model.var <- list() # grouped by a pack of AEQ subscales
for(AEQ in AEQs){ # reformulate(c("`A var`", "B"), "Y") AEQ ~ group + baseline.var
if(var.types[i] == "NA"){
formula <- as.formula(paste(AEQ,"s(group, bs = 'fs')", sep = " ~ "))
title <- sprintf("Response variable: %s\nPredictor Variable: group",AEQ)
}
else if(var.types[i] == "exam.score_study.time"){ # numeric
formula <- as.formula(
paste(AEQ,
paste("s(group, bs = 'fs')",paste0("s(",baseline.var,")"), sep = " + "), # use default setting
sep = " ~ "))
title <- sprintf("Response variable: %s\nPredictor Variables: group and %s\n",AEQ, baseline.var)
}
else{ # categorical: nominal and ordinal
formula <- as.formula(
paste(AEQ,
paste("s(group, bs = 'fs')",paste0("s(",baseline.var,", bs = 'fs')"), sep = " + "),
sep = " ~ "))
title <- sprintf("Response variable: %s\nPredictor Variables: group and %s\n",AEQ, baseline.var)
}
model <- mgcv::gam(formula, data = regression_unnest %>% na.omit(cols="AEQ"))
model.sum <- summary(model)
coefficient <- round(model$coefficients,2) # estimate, effect, slope of each value of smoothed variable
# p.coeff.smooth <- round(model$p.coeff,3) # p.value of each value of smoothed variable
sig.smooth.table <- model.sum$s.table # Approximate significance of smooth terms; combined effect of a variable (use F statistics, related to ANOVA)
# print(model.sum$p.table)
isSig <- ifelse(model.sum$s.table[,4] <= 0.05,T,F)
if(TRUE %in% isSig){
report <- append(report, paste0(strrep("*",times=60),"\n",title,"Approximate significance of smooth terms:\n"))
report <- append(report, paste0(capture.output(print(sig.smooth.table)),"\n"))
report <- append(report, "Coefficients:\n")
report <- append(report, paste0(capture.output(print(coefficient)),"\n"))
report <- append(report, "\n\n")
}
models.df[nrow(models.df) + 1,] = list(var.types[i], baseline.var, AEQ, list(coefficient), list(sig.smooth.table))
}
# names(model.group)[length(model.group)] <- baseline.var
}
i = i + 1
# names(models)[length(models)] <- baseline.var.type
}
report.regression <- report
cat(report.regression,sep = "")
## ************************************************************
## Response variable: AEQ.6.anxiety
## Predictor Variables: group and gender
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.963646e-11 3 4.080340e-12 0.536745919
## s(gender) 8.924445e-01 1 8.297529e+00 0.008117999
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(gender).1 s(gender).2
## 15.67 0.00 0.00 0.00 3.59 -3.59
##
##
## ************************************************************
## Response variable: AEQ.6.hopeless
## Predictor Variables: group and gender
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.537428e-11 3 1.575858e-12 0.72813180
## s(gender) 9.206092e-01 1 1.159592e+01 0.00291672
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(gender).1 s(gender).2
## 11.23 0.00 0.00 0.00 3.08 -3.08
##
##
## ************************************************************
## Response variable: AEQ.7.shame
## Predictor Variables: group and ethnicity
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.005383 2 1.275182 0.18833951
## s(ethnicity) 1.760722 2 6.137339 0.01274067
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(ethnicity).1
## 9.80 -0.45 -0.66 1.11 4.83
## s(ethnicity).2 s(ethnicity).3
## -1.67 -3.15
##
##
## ************************************************************
## Response variable: AEQ.6.anxiety
## Predictor Variables: group and meditation.times
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.7103534 2 0.7138021 0.28588164
## s(meditation.times) 1.7697671 2 5.4975622 0.01753378
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 18.03 -0.68 -0.51
## s(group).3 s(meditation.times).1 s(meditation.times).2
## 1.20 7.53 -2.61
## s(meditation.times).3
## -4.92
##
##
## ************************************************************
## Response variable: AEQ.14.pride
## Predictor Variables: group and meditation.times
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.125754e-10 3 3.219509e-11 0.43984921
## s(meditation.times) 1.233176e+00 2 2.478647e+00 0.04073988
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 12.99 0.00 0.00
## s(group).3 s(meditation.times).1 s(meditation.times).2
## 0.00 1.14 -2.08
## s(meditation.times).3
## 0.94
##
##
## ************************************************************
## Response variable: AEQ.14.anger
## Predictor Variables: group and concern.2
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.411756e-11 3 1.972596e-12 0.66472807
## s(concern.2) 1.655418e+00 2 3.657819e+00 0.03255441
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(concern.2).1
## 7.11 0.00 0.00 0.00 -2.05
## s(concern.2).2 s(concern.2).3
## 0.01 2.04
##
##
## ************************************************************
## Response variable: AEQ.6.hope
## Predictor Variables: group and concern.4
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 3.276738e-11 3 4.792431e-12 0.64941233
## s(concern.4) 1.514079e+00 2 3.218125e+00 0.03613567
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(concern.4).1
## 11.03 0.00 0.00 0.00 0.75
## s(concern.4).2 s(concern.4).3
## 2.94 -3.69
##
##
## ************************************************************
## Response variable: AEQ.14.anger
## Predictor Variables: group and concern.6
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.1649865 2 0.1127946 0.35435630
## s(concern.6) 1.7128591 2 4.1515520 0.02871288
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(concern.6).1
## 8.50 -0.10 0.11 -0.01 3.07
## s(concern.6).2 s(concern.6).3
## -0.33 -2.74
##
##
## ************************************************************
## Response variable: AEQ.6.hope
## Predictor Variables: group and religious.level
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.873650e-10 3 5.559132e-11 0.42878503
## s(religious.level) 1.156671e+00 3 1.355876e+00 0.04686367
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 11.73 0.00 0.00
## s(group).3 s(religious.level).1 s(religious.level).2
## 0.00 -2.73 0.92
## s(religious.level).3 s(religious.level).4
## 1.97 -0.16
##
##
## ************************************************************
## Response variable: AEQ.6.joy
## Predictor Variables: group and religious.level
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.7090335 2 0.4981796 0.29359816
## s(religious.level) 1.6078022 3 2.3868580 0.02545027
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 5.92 -0.36 0.03
## s(group).3 s(religious.level).1 s(religious.level).2
## 0.33 -1.15 -0.70
## s(religious.level).3 s(religious.level).4
## 1.34 0.52
##
##
## ************************************************************
## Response variable: AEQ.7.pride
## Predictor Variables: group and religious.level
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.775846 2 7.956171 0.01848977
## s(religious.level) 2.534293 3 7.525949 0.01462344
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 11.32 -0.53 2.55
## s(group).3 s(religious.level).1 s(religious.level).2
## -2.02 -1.11 -3.74
## s(religious.level).3 s(religious.level).4
## 2.38 2.46
##
##
## ************************************************************
## Response variable: AEQ.14.pride
## Predictor Variables: group and religious.level
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 3.541847e-11 3 8.767240e-12 0.49166972
## s(religious.level) 2.062603e+00 3 2.368111e+00 0.04729275
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 13.92 0.00 0.00
## s(group).3 s(religious.level).1 s(religious.level).2
## 0.00 -3.04 3.23
## s(religious.level).3 s(religious.level).4
## -0.76 0.57
##
##
## ************************************************************
## Response variable: AEQ.difference.joy
## Predictor Variables: group and religious.level
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.518268e-09 2 8.706893e-10 0.34431256
## s(religious.level) 2.196038e+00 3 3.393932e+00 0.01943028
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 1.55 0.00 0.00
## s(group).3 s(religious.level).1 s(religious.level).2
## 0.00 -1.72 -0.91
## s(religious.level).3 s(religious.level).4
## 0.60 2.03
##
##
## ************************************************************
## Response variable: AEQ.difference.pride
## Predictor Variables: group and religious.level
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 4.238554e-11 3 1.019828e-11 0.49939220
## s(religious.level) 1.318848e+00 3 1.636651e+00 0.03936183
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 0.98 0.00 0.00
## s(group).3 s(religious.level).1 s(religious.level).2
## 0.00 -1.65 0.53
## s(religious.level).3 s(religious.level).4
## 0.60 0.53
##
##
## ************************************************************
## Response variable: AEQ.difference.relief
## Predictor Variables: group and religious.level
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.1212395 2 0.07899692 0.31850355
## s(religious.level) 2.5952509 3 3.21635525 0.04057952
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 0.91 -0.03 0.05
## s(group).3 s(religious.level).1 s(religious.level).2
## -0.02 -1.07 3.16
## s(religious.level).3 s(religious.level).4
## -0.61 -1.47
##
##
## ************************************************************
## Response variable: AEQ.7.relief
## Predictor Variables: group and relax
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 2.517715e-11 3 3.146659e-12 0.69179466
## s(relax) 2.654209e+00 3 3.834007e+00 0.02537426
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(relax).1 s(relax).2
## 12.91 0.00 0.00 0.00 -7.18 2.99
## s(relax).3 s(relax).4
## 0.09 4.10
##
##
## ************************************************************
## Response variable: AEQ.6.hopeless
## Predictor Variables: group and improve.score
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.441772 2 2.850944 0.08649617
## s(improve.score) 1.555476 2 5.992362 0.01461130
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 12.59 -0.03 -2.17 2.20
## s(improve.score).1 s(improve.score).2 s(improve.score).3
## -2.74 -2.10 4.84
##
##
## ************************************************************
## Response variable: AEQ.difference.joy
## Predictor Variables: group and improve.score
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.0749976 2 0.0444193 0.33873609
## s(improve.score) 1.6180180 2 3.7230578 0.02996132
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 0.71 0.04 -0.03 -0.01
## s(improve.score).1 s(improve.score).2 s(improve.score).3
## -1.75 -0.46 2.21
##
##
## ************************************************************
## Response variable: AEQ.difference.joy
## Predictor Variables: group and race
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.929364 2 26.00661 0.01423183
## s(race) 1.694602 3 67.10021 0.02271690
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(race).1 s(race).2
## 0.17 2.03 -1.66 -0.38 1.36 0.47
## s(race).3 s(race).4
## -0.33 -1.50
##
##
## ************************************************************
## Response variable: AEQ.7.pride
## Predictor Variables: group and major
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.536758e+00 2 2.838733e+00 0.04536525
## s(major) 4.079446e-10 9 4.678239e-11 0.41872698
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(major).1 s(major).2
## 11.47 0.03 1.82 -1.84 0.00 0.00
## s(major).3 s(major).4 s(major).5 s(major).6 s(major).7 s(major).8
## 0.00 0.00 0.00 0.00 0.00 0.00
## s(major).9 s(major).10
## 0.00 0.00
##
##
## ************************************************************
## Response variable: AEQ.6.joy
## Predictor Variables: group and religion
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.6216313 2 0.4690953 0.29604160
## s(religion) 2.8145536 4 3.3254006 0.01670994
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(religion).1
## 5.70 -0.31 0.07 0.23 2.44
## s(religion).2 s(religion).3 s(religion).4 s(religion).5
## -2.28 0.37 -0.89 0.37
##
##
## ************************************************************
## Response variable: AEQ.7.pride
## Predictor Variables: group and religion
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.554670 2 3.457742 0.045157181
## s(religion) 1.809944 4 2.818659 0.008763601
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(religion).1
## 12.29 -0.06 1.65 -1.59 1.43
## s(religion).2 s(religion).3 s(religion).4 s(religion).5
## 0.69 -0.08 -2.29 0.26
##
##
## ************************************************************
## Response variable: AEQ.difference.anxiety
## Predictor Variables: group and religion
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.1250956 2 0.08206911 0.32255156
## s(religion) 3.3844217 4 3.63593183 0.02266672
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(religion).1
## -3.40 0.09 -0.07 -0.03 3.11
## s(religion).2 s(religion).3 s(religion).4 s(religion).5
## -6.70 -1.33 1.54 3.38
##
##
## ************************************************************
## Response variable: AEQ.difference.joy
## Predictor Variables: group and religion
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.107924 2 1.656441 0.105493556
## s(religion) 2.079555 4 3.294438 0.006621289
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(religion).1
## 1.47 0.65 -0.43 -0.22 1.02
## s(religion).2 s(religion).3 s(religion).4 s(religion).5
## 0.39 -0.45 -1.71 0.75
##
##
## ************************************************************
## Response variable: AEQ.difference.anxiety
## Predictor Variables: group and religion.family
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 2.656787e-11 3 6.852469e-12 0.47140680
## s(religion.family) 3.473816e+00 4 4.089946e+00 0.01295566
## Coefficients:
## (Intercept) s(group).1 s(group).2
## -3.94 0.00 0.00
## s(group).3 s(religion.family).1 s(religion.family).2
## 0.00 3.51 -6.27
## s(religion.family).3 s(religion.family).4 s(religion.family).5
## -0.83 1.45 2.13
##
##
## ************************************************************
## Response variable: AEQ.difference.shame
## Predictor Variables: group and religion.family
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 2.333945e-11 3 4.244912e-12 0.58225320
## s(religion.family) 3.345417e+00 4 2.830841e+00 0.03927289
## Coefficients:
## (Intercept) s(group).1 s(group).2
## -3.88 0.00 0.00
## s(group).3 s(religion.family).1 s(religion.family).2
## 0.00 3.27 -7.37
## s(religion.family).3 s(religion.family).4 s(religion.family).5
## -1.54 1.10 4.53
##
##
## ************************************************************
## Response variable: AEQ.6.anxiety
## Predictor Variables: group and meditation.type
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.147921 2 1.534574 0.14022336
## s(meditation.type) 2.992005 4 2.537193 0.04021984
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 17.55 -0.90 -1.04
## s(group).3 s(meditation.type).1 s(meditation.type).2
## 1.94 0.29 -2.14
## s(meditation.type).3 s(meditation.type).4 s(meditation.type).5
## -2.99 6.04 -1.20
##
##
## ************************************************************
## Response variable: AEQ.7.shame
## Predictor Variables: group and meditation.type
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.709489e+00 2 4.155214e+00 0.01905394
## s(meditation.type) 3.977563e-11 5 5.215105e-12 0.60223393
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 9.42 -1.03 -2.19
## s(group).3 s(meditation.type).1 s(meditation.type).2
## 3.22 0.00 0.00
## s(meditation.type).3 s(meditation.type).4 s(meditation.type).5
## 0.00 0.00 0.00
##
##
## ************************************************************
## Response variable: AEQ.14.anger
## Predictor Variables: group and meditation.type
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.9945681 2 1.403846 0.11275876
## s(meditation.type) 2.7329805 4 2.582344 0.02595243
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 7.60 -0.73 0.27
## s(group).3 s(meditation.type).1 s(meditation.type).2
## 0.46 0.70 -0.57
## s(meditation.type).3 s(meditation.type).4 s(meditation.type).5
## -0.95 -1.38 2.21
##
##
## ************************************************************
## Response variable: AEQ.14.shame
## Predictor Variables: group and meditation.type
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.625554 2 3.955241 0.06164407
## s(meditation.type) 3.298442 4 4.229381 0.01739870
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 8.22 -1.02 -0.98
## s(group).3 s(meditation.type).1 s(meditation.type).2
## 2.01 1.46 -0.17
## s(meditation.type).3 s(meditation.type).4 s(meditation.type).5
## -1.23 -3.50 3.44
##
##
## ************************************************************
## Response variable: AEQ.difference.hope
## Predictor Variables: group and meditation.type
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.570956e+00 2 5.132338e+00 0.006582657
## s(meditation.type) 1.008619e-10 5 1.690587e-11 0.492909946
## Coefficients:
## (Intercept) s(group).1 s(group).2
## -0.07 2.15 -1.18
## s(group).3 s(meditation.type).1 s(meditation.type).2
## -0.97 0.00 0.00
## s(meditation.type).3 s(meditation.type).4 s(meditation.type).5
## 0.00 0.00 0.00
##
##
## ************************************************************
## Response variable: AEQ.difference.joy
## Predictor Variables: group and meditation.type
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.355606e+00 2 2.828541e+00 0.02992754
## s(meditation.type) 8.357097e-11 5 9.980704e-12 0.64139394
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 0.55 1.11 -0.66
## s(group).3 s(meditation.type).1 s(meditation.type).2
## -0.45 0.00 0.00
## s(meditation.type).3 s(meditation.type).4 s(meditation.type).5
## 0.00 0.00 0.00
##
##
## ************************************************************
## Response variable: AEQ.6.anxiety
## Predictor Variables: group and credits
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.617952e-11 3 3.252343e-12 0.55178696
## s(credits) 5.145533e+00 6 3.817766e+00 0.01588944
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(credits).1 s(credits).2
## 17.97 0.00 0.00 0.00 -1.57 15.17
## s(credits).3 s(credits).4 s(credits).5 s(credits).6 s(credits).7
## 0.03 -5.14 -0.89 -3.64 -3.96
##
##
## ************************************************************
## Response variable: AEQ.6.hope
## Predictor Variables: group and credits
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.7607605 2 1.024785 0.40204462
## s(credits) 5.1511350 6 3.302533 0.04885178
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(credits).1 s(credits).2
## 11.76 -0.40 1.10 -0.70 -0.30 13.73
## s(credits).3 s(credits).4 s(credits).5 s(credits).6 s(credits).7
## -5.16 -0.34 -3.39 -3.03 -1.51
##
##
## ************************************************************
## Response variable: AEQ.6.hopeless
## Predictor Variables: group and credits
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.229261 2 1.696196 0.15505283
## s(credits) 2.236800 6 1.416187 0.03488754
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(credits).1 s(credits).2
## 10.87 0.43 -1.83 1.40 -1.13 -0.07
## s(credits).3 s(credits).4 s(credits).5 s(credits).6 s(credits).7
## 0.18 -1.42 -0.79 3.49 -0.27
##
##
## ************************************************************
## Response variable: AEQ.6.joy
## Predictor Variables: group and credits
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.607523e-11 3 5.773268e-14 0.98863051
## s(credits) 5.081982e+00 6 3.050695e+00 0.03182184
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(credits).1 s(credits).2
## 5.96 0.00 0.00 0.00 0.03 4.73
## s(credits).3 s(credits).4 s(credits).5 s(credits).6 s(credits).7
## -3.10 -0.66 -0.57 -1.24 0.81
##
##
## ************************************************************
## Response variable: AEQ.14.pride
## Predictor Variables: group and age
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.220416 2 1.475823 0.174047799
## s(age) 1.444133 2 6.381231 0.007180118
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(age).1 s(age).2
## 11.61 0.12 1.02 -1.14 -2.06 2.42
## s(age).3
## -0.36
##
##
## ************************************************************
## Response variable: AEQ.14.relief
## Predictor Variables: group and age
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.748432e-11 3 1.894520e-12 0.72515335
## s(age) 1.103376e+00 2 2.112201e+00 0.04666311
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(age).1 s(age).2
## 14.33 0.00 0.00 0.00 -1.54 2.02
## s(age).3
## -0.48
##
##
## ************************************************************
## Response variable: AEQ.7.pride
## Predictor Variables: group and GPA.previous
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 2.00002 2 2.190591e+02 0.01553894
## s(GPA.previous) 10.46619 12 8.420023e+11 0.02479328
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 11.10 -1.37 3.43 -2.06
## s(GPA.previous).1 s(GPA.previous).2 s(GPA.previous).3 s(GPA.previous).4
## -0.04 3.75 -1.52 2.69
## s(GPA.previous).5 s(GPA.previous).6 s(GPA.previous).7 s(GPA.previous).8
## 2.91 -4.01 1.29 -0.04
## s(GPA.previous).9 s(GPA.previous).10 s(GPA.previous).11 s(GPA.previous).12
## 3.75 -2.22 -1.06 -0.64
## s(GPA.previous).13
## -4.86
##
##
## ************************************************************
## Response variable: AEQ.7.shame
## Predictor Variables: group and GPA.previous
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.7393867 2 2.069981 0.22943482
## s(GPA.previous) 9.1638228 12 3.970002 0.03952446
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 8.49 0.15 -0.72 0.58
## s(GPA.previous).1 s(GPA.previous).2 s(GPA.previous).3 s(GPA.previous).4
## 6.65 -1.96 3.99 -0.60
## s(GPA.previous).5 s(GPA.previous).6 s(GPA.previous).7 s(GPA.previous).8
## 1.41 3.49 -2.80 -1.54
## s(GPA.previous).9 s(GPA.previous).10 s(GPA.previous).11 s(GPA.previous).12
## -1.96 0.18 -3.58 -3.45
## s(GPA.previous).13
## 0.18
##
##
## ************************************************************
## Response variable: AEQ.difference.anger
## Predictor Variables: group and GPA.previous
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 2.879492e-12 3 1.942965e-14 0.979456033
## s(GPA.previous) 1.139887e+01 12 1.787596e+01 0.002247874
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## -2.36 0.00 0.00 0.00
## s(GPA.previous).1 s(GPA.previous).2 s(GPA.previous).3 s(GPA.previous).4
## -10.01 -0.60 4.11 2.78
## s(GPA.previous).5 s(GPA.previous).6 s(GPA.previous).7 s(GPA.previous).8
## -5.95 -4.50 -6.25 4.11
## s(GPA.previous).9 s(GPA.previous).10 s(GPA.previous).11 s(GPA.previous).12
## 0.34 5.05 0.84 3.16
## s(GPA.previous).13
## 6.93
##
##
## ************************************************************
## Response variable: AEQ.difference.pride
## Predictor Variables: group and GPA.previous
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 2.185669e-11 3 5.553739e-12 0.49026982
## s(GPA.previous) 7.284372e+00 12 1.893067e+00 0.04810839
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 0.17 0.00 0.00 0.00
## s(GPA.previous).1 s(GPA.previous).2 s(GPA.previous).3 s(GPA.previous).4
## -0.09 1.58 -1.21 3.11
## s(GPA.previous).5 s(GPA.previous).6 s(GPA.previous).7 s(GPA.previous).8
## -0.48 -2.27 -0.09 -0.09
## s(GPA.previous).9 s(GPA.previous).10 s(GPA.previous).11 s(GPA.previous).12
## 2.14 -1.77 -1.20 -1.21
## s(GPA.previous).13
## 1.58
##
##
## ************************************************************
## Response variable: AEQ.6.anxiety
## Predictor Variables: group and grade.predict
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.430244 2 75.48078 0.090016847
## s(grade.predict) 8.548368 9 29.91289 0.004968139
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 16.29 1.39 -2.28 0.88
## s(grade.predict).1 s(grade.predict).2 s(grade.predict).3 s(grade.predict).4
## -1.88 0.32 5.66 -2.32
## s(grade.predict).5 s(grade.predict).6 s(grade.predict).7 s(grade.predict).8
## 1.60 18.74 -8.21 -0.96
## s(grade.predict).9 s(grade.predict).10
## -7.26 -5.69
##
##
## ************************************************************
## Response variable: AEQ.6.hope
## Predictor Variables: group and grade.predict
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.3064873 2 0.2416603 0.3962420
## s(grade.predict) 6.8992217 9 2.5714508 0.0498696
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 11.36 -0.33 0.34 -0.01
## s(grade.predict).1 s(grade.predict).2 s(grade.predict).3 s(grade.predict).4
## -2.65 -5.21 -5.48 -0.33
## s(grade.predict).5 s(grade.predict).6 s(grade.predict).7 s(grade.predict).8
## -0.72 11.86 -0.02 -0.50
## s(grade.predict).9 s(grade.predict).10
## -0.02 3.06
##
##
## ************************************************************
## Response variable: AEQ.6.joy
## Predictor Variables: group and grade.predict
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 7.464011e-11 3 2.162902e-11 0.43593760
## s(grade.predict) 7.154198e+00 9 2.842487e+00 0.03579758
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 5.75 0.00 0.00 0.00
## s(grade.predict).1 s(grade.predict).2 s(grade.predict).3 s(grade.predict).4
## -0.22 -1.87 -2.78 -0.69
## s(grade.predict).5 s(grade.predict).6 s(grade.predict).7 s(grade.predict).8
## 0.21 4.62 0.92 0.92
## s(grade.predict).9 s(grade.predict).10
## -1.30 0.18
##
##
## ************************************************************
## Response variable: AEQ.difference.hopeless
## Predictor Variables: group and grade.predict
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.3729402 2 0.4882811 0.33759603
## s(grade.predict) 7.5653496 9 3.9959974 0.02444592
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## -0.90 0.32 -0.19 -0.13
## s(grade.predict).1 s(grade.predict).2 s(grade.predict).3 s(grade.predict).4
## 1.60 -0.42 -7.95 -0.15
## s(grade.predict).5 s(grade.predict).6 s(grade.predict).7 s(grade.predict).8
## 3.42 -3.99 -1.14 3.29
## s(grade.predict).9 s(grade.predict).10
## 3.67 1.68
##
##
## ************************************************************
## Response variable: AEQ.6.anxiety
## Predictor Variables: group and productivity6
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 8.774577e-10 3 3.207978e-10 0.3615621055
## s(productivity6) 4.785494e+00 5 1.035719e+01 0.0005892683
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 17.06 0.00 0.00 0.00
## s(productivity6).1 s(productivity6).2 s(productivity6).3 s(productivity6).4
## -7.53 -0.05 0.44 18.64
## s(productivity6).5 s(productivity6).6
## -3.03 -8.46
##
##
## ************************************************************
## Response variable: AEQ.6.hope
## Predictor Variables: group and productivity6
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.300778 2 7.172321 0.137565173
## s(productivity6) 4.653697 5 8.370785 0.006635131
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 12.90 -0.01 1.66 -1.65
## s(productivity6).1 s(productivity6).2 s(productivity6).3 s(productivity6).4
## -1.70 -6.24 -5.33 15.15
## s(productivity6).5 s(productivity6).6
## -3.19 1.31
##
##
## ************************************************************
## Response variable: AEQ.6.joy
## Predictor Variables: group and productivity6
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.968938e-11 3 2.702484e-12 0.66371958
## s(productivity6) 4.384781e+00 5 4.173470e+00 0.01259658
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 5.96 0.00 0.00 0.00
## s(productivity6).1 s(productivity6).2 s(productivity6).3 s(productivity6).4
## 0.85 -2.22 -1.62 4.93
## s(productivity6).5 s(productivity6).6
## -0.33 -1.60
##
##
## ************************************************************
## Response variable: AEQ.14.anger
## Predictor Variables: group and productivity6
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 0.6524268 2 0.5417918 0.25320991
## s(productivity6) 2.1429294 5 1.5496847 0.03465215
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 6.43 -0.49 0.27 0.22
## s(productivity6).1 s(productivity6).2 s(productivity6).3 s(productivity6).4
## -0.27 -1.04 0.11 -0.19
## s(productivity6).5 s(productivity6).6
## 2.18 -0.79
##
##
## ************************************************************
## Response variable: AEQ.14.pride
## Predictor Variables: group and productivity6
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.967035 2 53.58624 0.002265159
## s(productivity6) 4.359018 5 85.43736 0.004438901
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 13.71 1.43 2.26 -3.69
## s(productivity6).1 s(productivity6).2 s(productivity6).3 s(productivity6).4
## -0.11 0.39 -5.53 4.99
## s(productivity6).5 s(productivity6).6
## -3.09 3.36
##
##
## ************************************************************
## Response variable: AEQ.6.hopeless
## Predictor Variables: group and productivity7
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 9.758472e-12 3 1.105974e-12 0.71331338
## s(productivity7) 2.502778e+00 3 2.872017e+00 0.04916358
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 12.30 0.00 0.00 0.00
## s(productivity7).1 s(productivity7).2 s(productivity7).3 s(productivity7).4
## -4.88 4.97 -1.19 1.10
##
##
## ************************************************************
## Response variable: AEQ.7.anger
## Predictor Variables: group and productivity7
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 5.245249e-11 3 1.384312e-11 0.46802089
## s(productivity7) 2.793599e+00 3 4.790408e+00 0.01461879
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3
## 8.41 0.00 0.00 0.00
## s(productivity7).1 s(productivity7).2 s(productivity7).3 s(productivity7).4
## -0.15 3.85 -1.53 -2.17
##
##
## ************************************************************
## Response variable: AEQ.7.anger
## Predictor Variables: group and exam.score
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 2.629128e-11 3.00000 2.534469e-12 0.73127260
## s(exam.score) 2.076956e+00 2.58593 4.043894e+00 0.03767037
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(exam.score).1
## 7.29 0.00 0.00 0.00 -0.41
## s(exam.score).2 s(exam.score).3 s(exam.score).4 s(exam.score).5 s(exam.score).6
## 0.45 -0.13 0.52 -0.20 0.50
## s(exam.score).7 s(exam.score).8 s(exam.score).9
## 0.00 -3.02 0.95
##
##
## ************************************************************
## Response variable: AEQ.difference.anger
## Predictor Variables: group and exam.score
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.232233 2 1.827157 0.084401806
## s(exam.score) 1.000000 1 10.773325 0.005454788
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(exam.score).1
## -3.07 2.01 -0.45 -1.56 0.00
## s(exam.score).2 s(exam.score).3 s(exam.score).4 s(exam.score).5 s(exam.score).6
## 0.00 0.00 0.00 0.00 0.00
## s(exam.score).7 s(exam.score).8 s(exam.score).9
## 0.00 0.00 3.37
##
##
## ************************************************************
## Response variable: AEQ.difference.relief
## Predictor Variables: group and exam.score
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.308655 2.000000 3.462043 0.04608096
## s(exam.score) 8.427342 8.897769 4.165225 0.05169583
## Coefficients:
## (Intercept) s(group).1 s(group).2 s(group).3 s(exam.score).1
## 0.21 -0.86 0.99 -0.13 -7.33
## s(exam.score).2 s(exam.score).3 s(exam.score).4 s(exam.score).5 s(exam.score).6
## 4.97 1.09 2.24 2.11 0.05
## s(exam.score).7 s(exam.score).8 s(exam.score).9
## -1.41 -4.74 11.03
##
##
## ************************************************************
## Response variable: AEQ.7.pride
## Predictor Variables: group and study.time.total
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.700853 2.000000 4.727285 0.01715818
## s(study.time.total) 1.276414 1.494543 5.235040 0.02018309
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 10.98 0.39 2.05
## s(group).3 s(study.time.total).1 s(study.time.total).2
## -2.43 -0.03 0.21
## s(study.time.total).3 s(study.time.total).4 s(study.time.total).5
## 0.06 0.15 0.06
## s(study.time.total).6 s(study.time.total).7 s(study.time.total).8
## -0.17 -0.12 0.77
## s(study.time.total).9
## 1.72
##
##
## ************************************************************
## Response variable: AEQ.7.shame
## Predictor Variables: group and study.time.total
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.690449 2 3.970114 0.02918300
## s(study.time.total) 1.000000 1 4.388741 0.05633225
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 8.93 -1.13 -2.17
## s(group).3 s(study.time.total).1 s(study.time.total).2
## 3.30 0.00 0.00
## s(study.time.total).3 s(study.time.total).4 s(study.time.total).5
## 0.00 0.00 0.00
## s(study.time.total).6 s(study.time.total).7 s(study.time.total).8
## 0.00 0.00 0.00
## s(study.time.total).9
## -1.75
##
##
## ************************************************************
## Response variable: AEQ.14.anger
## Predictor Variables: group and study.time.total
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.893573 2.000000 5.157426 0.04035853
## s(study.time.total) 8.035068 8.737403 4.698256 0.03564614
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 7.66 -2.26 -0.24
## s(group).3 s(study.time.total).1 s(study.time.total).2
## 2.50 11.81 -42.79
## s(study.time.total).3 s(study.time.total).4 s(study.time.total).5
## 1.50 -26.48 -3.28
## s(study.time.total).6 s(study.time.total).7 s(study.time.total).8
## 24.32 14.47 -61.53
## s(study.time.total).9
## 20.60
##
##
## ************************************************************
## Response variable: AEQ.14.shame
## Predictor Variables: group and study.time.total
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.979829 2.000000 4.001508 0.04857260
## s(study.time.total) 6.214019 7.304922 2.748287 0.08876767
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 8.09 -2.51 -1.51
## s(group).3 s(study.time.total).1 s(study.time.total).2
## 4.03 12.89 -36.87
## s(study.time.total).3 s(study.time.total).4 s(study.time.total).5
## 6.73 -22.62 -4.20
## s(study.time.total).6 s(study.time.total).7 s(study.time.total).8
## 15.71 12.09 -51.19
## s(study.time.total).9
## 23.61
##
##
## ************************************************************
## Response variable: AEQ.difference.hope
## Predictor Variables: group and study.time.total
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.486005 2 3.011461 0.04093086
## s(study.time.total) 1.000000 1 3.736281 0.07373060
## Coefficients:
## (Intercept) s(group).1 s(group).2
## -0.46 1.85 -0.41
## s(group).3 s(study.time.total).1 s(study.time.total).2
## -1.44 0.00 0.00
## s(study.time.total).3 s(study.time.total).4 s(study.time.total).5
## 0.00 0.00 0.00
## s(study.time.total).6 s(study.time.total).7 s(study.time.total).8
## 0.00 0.00 0.00
## s(study.time.total).9
## 1.26
##
##
## ************************************************************
## Response variable: AEQ.difference.joy
## Predictor Variables: group and study.time.total
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(group) 1.380504 2 2.477424 0.055260845
## s(study.time.total) 1.000000 1 10.330455 0.006246406
## Coefficients:
## (Intercept) s(group).1 s(group).2
## 0.51 1.03 -0.35
## s(group).3 s(study.time.total).1 s(study.time.total).2
## -0.68 0.00 0.00
## s(study.time.total).3 s(study.time.total).4 s(study.time.total).5
## 0.00 0.00 0.00
## s(study.time.total).6 s(study.time.total).7 s(study.time.total).8
## 0.00 0.00 0.00
## s(study.time.total).9
## 1.30
report <- NULL
datatable(
models.df, filter = 'top', extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)