The Importance of Client–Canine Contact in Canine-Assisted Interventions: A Randomized Controlled Trial

This document provides the R code and output for all analyses, as well as code for reproducing tables and figures in the manuscript. In addition to allowing for reviewing of the code, it also provides additional information that would not fit in the manuscript, such as figures and analyses conducted to evaluate test assumptions.

Analyses were conducted in R using Jupyter Hub.

Setup

Load Packages

library(tidyr)
library(purrr)
library(numform)
library(psych)
library(dplyr)
library(ggplot2)
library(ggtext)
library(forcats)
library(lmtest)
library(sandwich)

Attaching package: ‘dplyr’


The following object is masked from ‘package:numform’:

    collapse


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


Warning message:
“package ‘ggplot2’ was built under R version 4.0.3”

Attaching package: ‘ggplot2’


The following objects are masked from ‘package:psych’:

    %+%, alpha



Attaching package: ‘forcats’


The following object is masked from ‘package:numform’:

    as_factor


Loading required package: zoo


Attaching package: ‘zoo’


The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

User-defined Functions

format_num

This function has no impact on the analyses. I used it to format the numbers before putting them in a table.

format_num <- function(x) {
    x <- sprintf("%0.2f", x)
    gsub("-", "−", x)
}

t_test

This function does a paired samples t test and returns the result as a data.frame. Compares t2 to t1 and performs a directional test for the well-being measures and a two-sided test for the other measures.

t_test <- function(x, h = c("post", "follow-up")) {
    h <- match.arg(h)
    wb_vars <- c("FS", "PANAS-PA", "SCS", "SHS", "Composite")
    ib_vars <- c("Stress", "Homesick", "Lonely", "PANAS-NA")
    
    diff <- x$t2 - x$t1
    
    if (unique(x$measure) %in% wb_vars) {
        alt <- "greater"
    } else if (unique(x$measure) %in% ib_vars) {
        alt <- "less"
    } else {
        alt <- "two.sided"
    }
    
    result <- t.test(x = diff, alternative = alt)
    d <- psych::t2d(result$statistic, n1 = result$parameter + 1)
    d90 <- psych::d.ci(d, n1 = result$parameter + 1, alpha = .10)
    d95 <- psych::d.ci(d, n1 = result$parameter + 1)
    data.frame(
        t = abs(result$statistic),
        df = result$parameter,
        p = result$p.value,
        d = d95[[2]],
        d90.LL = d90[[1]],
        d90.UL = d90[[3]],
        d95.LL = d95[[1]],
        d95.UL = d95[[3]]
    )
}

ParallelAnalysis

This function conducts parallel analysis, which is generally considered the gold standard for determining the number of factors to retain in exploratory factor analysis. Used in the exploratory analyses.

ParallelAnalysis <- function(n_obs, n_vars, percentiles = c(.5, .95),
                             n_iterations = 1000) {
    rdta <- replicate(
        n = n_iterations,
        expr = eigen(
            x = cor(matrix(rnorm(n_obs * n_vars), ncol = n_vars)),
            symmetric = TRUE,
            only.values = TRUE)$values,
        simplify = "array"
    )

    pa_eigvals <- apply(rdta,
                        MARGIN = 1,
                        FUN = quantile,
                        probs = percentiles)

    return(t(pa_eigvals))
}

ctt and multi_ctt

These stand for Classical Test Theory. They just give pretty output of classical test theory reliability analysis.

ctt <- function(x, varnames, keys, drop = NULL, pretty_varnames = NULL) {
    if (!is.null(drop)) {
        include <- varnames[-which(varnames %in% drop)]
    } else {
        include <- varnames
    }
    
    if (is.null(pretty_varnames)) pretty_varnames <- varnames
    
    result <- suppressMessages(psych::alpha(x[, include], keys = keys))
    
    dt <- data.frame(
        variable = include,
        rdrop = f_num(result$item.stats$r.drop, digits = 2)
    )
    
    out <- data.frame(
        variable = varnames,
        Item = pretty_varnames
    )
    
    out <- left_join(
        x = out,
        y = dt,
        by = "variable"
    )
    
    list(rdrop = out[, c("Item", "rdrop")], alpha = f_num(result$total$raw_alpha, digits = 2))
}

multi_ctt <- function(x, varnames, keys, drop, pretty_varnames = NULL) {
    # Step 1
    out <- ctt(x, varnames, keys, drop = NULL, pretty_varnames)
    for (i in seq_along(drop)) {
        res <- ctt(x, varnames, keys, drop = drop[1:i], pretty_varnames)
        out[[1]] <- left_join(
            x = out[[1]],
            y = res[[1]],
            by = "Item"
        )
        out[[2]] <- c(out[[2]], res[[2]])
    }
    names(out[[1]]) <- c("Item", paste0("Step", seq_len(length(drop) + 1)))
    out[[1]] <- sapply(out[[1]], function(x) ifelse(is.na(x), "-", x))
    out
}

Data Import and Wrangling

The dataset being imported includes composite scores for all scales measures at pre-intervention, post-intervention, and at two-week follow-up.

dta <- read.csv("canine-contact-scored.csv")
dta$group <- factor(dta$group,
    levels = c("Handler", "Indirect", "Direct")
)
dta$time <- factor(
    dta$time,
    levels = c("t1", "t2")
)
head(dta)

A data.frame: 6 × 13

	pid	group	time	pets	stress	homesickness	loneliness	panas_pa	panas_na	shs	scs	flourishing	integration
	<int>	<fct>	<fct>	<chr>	<int>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<int>
1	1	Handler	t1	No	2	3	2.25	3.2	1.0	2.333333	3.900000	6.125	1
2	1	Handler	t2	No	2	3	1.70	3.8	1.2	3.333333	3.800000	6.000	2
3	2	Direct	t1	No	2	4	1.90	3.0	1.8	3.333333	5.150000	5.250	4
4	2	Direct	t2	No	1	2	1.60	3.2	1.0	4.000000	5.263158	6.000	4
5	3	Indirect	t1	No	4	3	2.25	2.8	1.6	2.666667	4.100000	5.375	3
6	3	Indirect	t2	No	3	2	2.25	3.0	1.6	3.333333	4.150000	5.375	3

Reshape the Data

This is just rearranging the data frame so that each measure is a column. This simplifies the code later by allowing us to apply a function to the data at each level of the measure column. The analyses are conducted using the dta_diff.

dta_longer <- pivot_longer(
    dta,
    cols = c(-pid, -group, -time, -pets),
    names_to = "measure",
    values_to = "score"
)

head(dta_longer)

dta_longer$measure <- factor(
    dta_longer$measure,
levels = c(
    "flourishing",
    "panas_pa",
    "scs",
    "shs",
    "integration",
    "stress",
    "homesickness",
    "loneliness",
    "panas_na"
  ),
  labels = c(
    "FS",
    "PANAS-PA",
    "SCS-R",
    "SHS",
    "Integration",
    "Stress",
    "Homesick",
    "Lonely",
    "PANAS-NA"
  )
)

dta_diff <- pivot_wider(
    dta_longer,
    names_from = "time",
    values_from = "score"
)

dta_diff$group.measure <- paste(dta_diff$group, dta_diff$measure, sep = ".")

dta_diff$diff21 <- dta_diff$t2 - dta_diff$t1

head(dta_diff)

A tibble: 6 × 6

pid	group	time	pets	measure	score
<int>	<fct>	<fct>	<chr>	<chr>	<dbl>
1	Handler	t1	No	stress	2.000000
1	Handler	t1	No	homesickness	3.000000
1	Handler	t1	No	loneliness	2.250000
1	Handler	t1	No	panas_pa	3.200000
1	Handler	t1	No	panas_na	1.000000
1	Handler	t1	No	shs	2.333333

A tibble: 6 × 8

pid	group	pets	measure	t1	t2	group.measure	diff21
<int>	<fct>	<chr>	<fct>	<dbl>	<dbl>	<chr>	<dbl>
1	Handler	No	Stress	2.000000	2.000000	Handler.Stress	0.00
1	Handler	No	Homesick	3.000000	3.000000	Handler.Homesick	0.00
1	Handler	No	Lonely	2.250000	1.700000	Handler.Lonely	-0.55
1	Handler	No	PANAS-PA	3.200000	3.800000	Handler.PANAS-PA	0.60
1	Handler	No	PANAS-NA	1.000000	1.200000	Handler.PANAS-NA	0.20
1	Handler	No	SHS	2.333333	3.333333	Handler.SHS	1.00

Scale Descriptives

These are descriptive values for each scale.

Note that the analyses were conducted on either change scores (H1 and H2) or the residual post-intervention score, after controlling for pre-test score. As such, the tests do not make any assumptions about the distributions of these measures. Values such as skew and kurtosis are provided for descriptive purposes.

round(describe(dta[dta$time == "t1", -1:-4]), 3)
round(describe(dta[dta$time == "t2", -1:-4]), 3)

A psych: 9 × 13

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
stress	1	284	3.208	0.971	3.000	3.206	1.483	1.00	5.0	4.00	-0.145	-0.619	0.058
homesickness	2	284	2.088	1.205	2.000	1.943	1.483	1.00	5.0	4.00	0.738	-0.652	0.071
loneliness	3	284	2.040	0.530	1.950	2.008	0.519	1.05	3.6	2.55	0.513	-0.246	0.031
panas_pa	4	284	3.289	0.683	3.400	3.312	0.593	1.20	5.0	3.80	-0.355	0.356	0.041
panas_na	5	284	1.979	0.705	1.800	1.911	0.593	1.00	4.0	3.00	0.710	-0.331	0.042
shs	6	284	3.311	0.876	3.333	3.348	0.988	1.00	5.0	4.00	-0.328	-0.334	0.052
scs	7	284	4.416	0.824	4.550	4.478	0.890	1.60	5.9	4.30	-0.735	0.354	0.049
flourishing	8	284	5.678	0.782	5.875	5.745	0.741	2.75	7.0	4.25	-0.877	0.780	0.046
integration	9	284	3.222	1.048	3.000	3.224	1.483	1.00	5.0	4.00	-0.229	-0.589	0.062

A psych: 9 × 13

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
stress	1	284	2.373	1.030	2.000	2.311	1.483	1.000	5.00	4.000	0.466	-0.453	0.061
homesickness	2	284	1.894	1.113	1.000	1.711	0.000	1.000	5.00	4.000	1.065	0.141	0.066
loneliness	3	284	1.914	0.535	1.850	1.889	0.593	1.000	3.55	2.550	0.400	-0.540	0.032
panas_pa	4	284	3.296	0.787	3.400	3.298	0.890	1.200	6.80	5.600	0.136	0.696	0.047
panas_na	5	284	1.613	0.671	1.400	1.500	0.593	1.000	4.40	3.400	1.418	1.837	0.040
shs	6	284	3.534	0.823	3.667	3.558	0.988	1.333	5.00	3.667	-0.285	-0.460	0.049
scs	7	284	4.534	0.867	4.700	4.601	0.890	0.000	6.00	6.000	-0.954	1.843	0.051
flourishing	8	284	5.742	0.780	5.875	5.803	0.741	2.750	7.00	4.250	-0.829	0.654	0.046
integration	9	283	3.353	1.026	3.000	3.374	1.483	1.000	5.00	4.000	-0.331	-0.407	0.061

Test Hypothesis 1

Research question: Does the intervention improve subjective well-being?

Hypothesis: Participants in each of the three conditions will report significant improvements in subjective well-being from pre-to-post intervention.

This was tested by conducting a paired samples t test comparing pre-to-post-intervention changes on each of the outcome measures. Directional tests were used for measures of well-being (flourishing, positive affect, social connectedness, happiness, and integration with campus community), and ill-being (stress, homesickness, loneliness, negative affect).

Descriptive Statistics

The paired samples t test is just a one-sample t test conducted on the change scores. So, the assumptions of the test are about the change scores. Below are descriptive statistics for the change scores.

desc_h1 <- pivot_wider(
    dta_diff[, c("pid", "group", "measure", "diff21")],
    names_from = measure,
    values_from = diff21
)

desc_h1 <- describe(desc_h1[, -1:-2])

round(desc_h1, 3)

A psych: 9 × 13

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
Stress	1	284	-0.835	0.953	-1.0	-0.842	1.483	-4.00	2.000	6.000	-0.065	0.397	0.057
Homesick	2	284	-0.194	0.854	0.0	-0.171	0.000	-3.00	2.000	5.000	-0.436	1.864	0.051
Lonely	3	284	-0.126	0.229	-0.1	-0.121	0.222	-1.15	0.650	1.800	-0.366	1.758	0.014
PANAS-PA	4	284	0.007	0.511	0.0	0.013	0.297	-1.80	2.600	4.400	0.089	3.167	0.030
PANAS-NA	5	284	-0.366	0.510	-0.2	-0.331	0.297	-3.00	1.200	4.200	-0.965	3.058	0.030
SHS	6	284	0.223	0.410	0.0	0.211	0.494	-1.00	1.667	2.667	0.397	0.712	0.024
SCS-R	7	284	0.118	0.412	0.1	0.124	0.232	-2.90	2.200	5.100	-1.490	18.206	0.024
FS	8	284	0.064	0.351	0.0	0.057	0.185	-1.00	1.375	2.375	0.351	1.211	0.021
Integration	9	283	0.131	0.695	0.0	0.154	0.000	-3.00	2.000	5.000	-0.432	3.614	0.041

Histograms

Histograms of the difference scores.

ggplot(dta_diff, aes(x = t2 - t1)) +
    facet_grid(rows = vars(measure), cols = vars(group)) +
    geom_histogram(bins = 15) +
    theme_minimal() +
    theme(strip.text.y = element_text(angle = 0))

Warning message:
“Removed 1 rows containing non-finite values (stat_bin).”

results1 <- lapply(split(dta_diff, dta_diff$group.measure), t_test, h = "post")

results1 <- bind_rows(results1, .id = "Group.Measure")

results1 <- separate(
    results1,
    Group.Measure,
    into = c("Group", "Measure"),
    sep = "\\."
)

results1$Group <- factor(results1$Group, levels = c("Handler", "Indirect", "Direct"))

results1$Measure <- factor(results1$Measure, levels = levels(dta_diff$measure))

results1 <- results1[order(results1$Group, results1$Measure), ]

rownames(results1) <- NULL

results1

A data.frame: 27 × 10

Group	Measure	t	df	p	d	d90.LL	d90.UL	d95.LL	d95.UL
<fct>	<fct>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
Handler	FS	0.2182676	93	4.138494e-01	0.02251258	-0.147221060	0.192128500	-0.17972863	0.22463347
Handler	PANAS-PA	0.9237790	93	8.210040e-01	-0.09528053	-0.265066876	0.075013623	-0.29764294	0.10759107
Handler	SCS-R	0.6288815	93	5.309683e-01	0.06486418	-0.105143219	0.234523801	-0.13767854	0.26705991
Handler	SHS	3.0553471	93	1.466669e-03	0.31513499	0.140464314	0.488160012	0.10717531	0.52147980
Handler	Integration	1.2642567	93	2.092973e-01	-0.13039812	-0.300428930	0.040329391	-0.33307007	0.07296862
Handler	Stress	5.4237260	93	2.299012e-07	-0.55941463	-0.740573555	-0.375483872	-0.77562323	-0.34059398
Handler	Homesick	1.7144738	93	4.488626e-02	-0.17683447	-0.347349252	-0.005378766	-0.38010840	0.02737496
Handler	Lonely	2.7803816	93	3.285873e-03	-0.28677446	-0.459156314	-0.112888118	-0.49233575	-0.07973336
Handler	PANAS-NA	5.2105989	93	5.639591e-07	-0.53743224	-0.717678843	-0.354509113	-0.75253489	-0.31979456
Indirect	FS	1.7564329	94	4.113561e-02	0.18020622	0.009598236	0.349865768	-0.02299259	0.38246182
Indirect	PANAS-PA	0.8215723	94	7.933004e-01	-0.08429154	-0.253129871	0.084992547	-0.28551930	0.11737819
Indirect	SCS-R	2.3151966	94	2.277800e-02	0.23753416	0.065772907	0.408054398	0.03299290	0.44084744
Indirect	SHS	5.0246541	94	1.198886e-06	0.51551863	0.334535463	0.693946801	0.30017196	0.72843575
Indirect	Integration	2.1619512	94	3.316112e-02	0.22181152	0.050395422	0.392068371	0.01766779	0.42480207
Indirect	Stress	7.1173652	94	1.081636e-10	-0.73022626	-0.918607898	-0.538436211	-0.95517926	-0.50217794
Indirect	Homesick	0.6279294	94	2.657864e-01	-0.06442419	-0.233188170	0.104681536	-0.26555211	0.13704467
Indirect	Lonely	4.1009251	94	4.371153e-05	-0.42074604	-0.595806703	-0.243553316	-0.62958244	-0.20984925
Indirect	PANAS-NA	5.9112079	94	2.714523e-08	-0.60647714	-0.788745574	-0.421274311	-0.82404005	-0.38617501
Direct	FS	3.3857567	94	5.186909e-04	0.34737131	0.172665993	0.520289965	0.13938803	0.55361121
Direct	PANAS-PA	2.3132034	94	1.144605e-02	0.23732966	0.065573017	0.407846322	0.03279376	0.44063859
Direct	SCS-R	5.5406018	94	2.736870e-07	0.56845375	0.385080859	0.749046699	0.35030115	0.78399046
Direct	SHS	7.8764449	94	2.897339e-12	0.80810620	0.611707895	1.000821510	0.57463984	1.03829348
Direct	Integration	5.0025133	93	2.662956e-06	0.51596985	0.333999147	0.695356229	0.29945002	0.73003238
Direct	Stress	15.4004846	94	9.510848e-28	-1.58005638	-1.830792184	-1.323604138	-1.88025161	-1.27592883
Direct	Homesick	4.0112455	94	6.057669e-05	-0.41154510	-0.586314895	-0.234690154	-0.62003108	-0.20103957
Direct	Lonely	9.7137481	94	3.730558e-16	-0.99660953	-1.201192067	-0.787744131	-1.24112917	-0.74848120
Direct	PANAS-NA	10.1735671	94	3.932684e-17	-1.04378596	-1.251615268	-0.831537743	-1.29222532	-0.79168106

Create a Pretty Table of the Results for H1

The code below simply formats the results so that they look nice and can be copied and pasted into MS Word to create a table, without having to manually enter the table data (which is time consuming and prone to human error).

tbl1 <- results1

tbl1$p <- f_num(tbl1$p, digits = 3)
tbl1$p <- ifelse(tbl1$p == ".000", "<.001", tbl1$p)

tbl1$df <- as.character(tbl1$df)

tbl1[, sapply(tbl1, is.numeric)] <- sapply(
    tbl1[, sapply(tbl1, is.numeric)],
    format_num
)

tbl1$Measure <- as.character(tbl1$Measure)

tbl1 <- split(tbl1[-1], tbl1$Group)

keepers1 <- c("Measure", "t", "df", "p", "d", "d95.LL", "d95.UL")

tbl1 <- rbind(
    c(names(tbl1)[[1]], rep("", length(keepers1) - 1)),
    tbl1[[1]][, keepers1],
    c(names(tbl1)[[2]], rep("", length(keepers1) - 1)),
    tbl1[[2]][, keepers1],
    c(names(tbl1)[[3]], rep("", length(keepers1) - 1)),
    tbl1[[3]][, keepers1]
)

rownames(tbl1) <- NULL

Copy and paste this output into MS Word, and then use Insert > Table > Convert Text to Table. There is still some formatting that needs to be done, but not much!

write.csv(tbl1, row.names = FALSE, quote = FALSE)

Measure,t,df,p,d,d95.LL,d95.UL
Handler,,,,,,
FS,0.22,93,.414,0.02,−0.18,0.22
PANAS-PA,0.92,93,.821,−0.10,−0.30,0.11
SCS-R,0.63,93,.531,0.06,−0.14,0.27
SHS,3.06,93,.001,0.32,0.11,0.52
Integration,1.26,93,.209,−0.13,−0.33,0.07
Stress,5.42,93,<.001,−0.56,−0.78,−0.34
Homesick,1.71,93,.045,−0.18,−0.38,0.03
Lonely,2.78,93,.003,−0.29,−0.49,−0.08
PANAS-NA,5.21,93,<.001,−0.54,−0.75,−0.32
Indirect,,,,,,
FS,1.76,94,.041,0.18,−0.02,0.38
PANAS-PA,0.82,94,.793,−0.08,−0.29,0.12
SCS-R,2.32,94,.023,0.24,0.03,0.44
SHS,5.02,94,<.001,0.52,0.30,0.73
Integration,2.16,94,.033,0.22,0.02,0.42
Stress,7.12,94,<.001,−0.73,−0.96,−0.50
Homesick,0.63,94,.266,−0.06,−0.27,0.14
Lonely,4.10,94,<.001,−0.42,−0.63,−0.21
PANAS-NA,5.91,94,<.001,−0.61,−0.82,−0.39
Direct,,,,,,
FS,3.39,94,.001,0.35,0.14,0.55
PANAS-PA,2.31,94,.011,0.24,0.03,0.44
SCS-R,5.54,94,<.001,0.57,0.35,0.78
SHS,7.88,94,<.001,0.81,0.57,1.04
Integration,5.00,93,<.001,0.52,0.30,0.73
Stress,15.40,94,<.001,−1.58,−1.88,−1.28
Homesick,4.01,94,<.001,−0.41,−0.62,−0.20
Lonely,9.71,94,<.001,−1.00,−1.24,−0.75
PANAS-NA,10.17,94,<.001,−1.04,−1.29,−0.79

Figure for H1

The code below creates the figure used in the manuscript (although using a different font).

ggplot(results1, aes(x = d, y = fct_rev(Measure))) +
    facet_wrap(~Group) +
    geom_vline(xintercept = 0, size = 1.0, colour = "#e9e9e9") +
    geom_linerange(aes(xmin = d95.LL, xmax = d95.UL)) +
    geom_linerange(aes(xmin = d90.LL, xmax = d90.UL), size = 1) +
    geom_point() +
    theme_minimal(base_family = "Fira Sans") +
    labs(
        x = "Cohen's *d<sub>z</sub>*",
        y = NULL
    ) +
    theme(
        axis.title.x = element_markdown(),
        strip.text = element_text(size = 12, hjust = 0, face = "bold")
    )

ggsave("h1.png", width = 6.5, height = 4, dpi = 320)

ANCOVAS

contrasts(dta_diff$group) <- cbind(
    handler_v_dog = c(-1, 0.5, 0.5),
    indirect_v_direct = c(0, -1, 1)
)

contrasts(dta_diff$group)

A matrix: 3 × 2 of type dbl

	handler_v_dog	indirect_v_direct
Handler	-1.0	0
Indirect	0.5	-1
Direct	0.5	1

models <- lapply(
    split(dta_diff, dta_diff$measure),
    function(x) aov(t2 ~ t1 + group, data = x)
)

corrected_models <- lapply(
    models,
    function(x) lmtest::coeftest(x, vcov. = sandwich::vcovHC(x))
)

Q–Q Plots

These plots show the relationship between the theoretical quantiles (i.e., normally distributed residuals) and the sample quantiles. If the sample residuals are normally distributed, the theoretical and sample quantiles should fall along the line.

plot(models$FS, which = 2, main = "FS")
plot(models$`PANAS-PA`, which = 2, main = "PANAS-PA")
plot(models$`SCS-R`, which = 2, main = "SCS-R")
plot(models$SHS, which = 2, main = "SHS")
plot(models$Integration, which = 2, main = "Integration")
plot(models$Stress, which = 2, main  = "Stress")
plot(models$Homesick, which = 2, main  = "Homesick")
plot(models$Lonely, which = 2, main  = "Lonely")
plot(models$`PANAS-NA`, which = 2, main  = "PANAS-NA")

results3 <- map2(
    models, corrected_models,
    ~ tibble(
        contrast = factor(1:2,
            levels = 1:2,
            labels = c("Handler v. Indirect and Direct", "Indirect v. Direct")
        ),
        n = length(.x$residuals),
        B = .y[3:4, "Estimate"],
        SE = .y[3:4, "Std. Error"],
        t = abs(.y[3:4, "t value"]),
        df = n - 4,
        p = .y[3:4, "Pr(>|t|)"],
        d = t2d(t = .y[3:4, "t value"], n = df),
        LL90 = d.ci(d = d, n = df, alpha = .10)[1:2, 1],
        UL90 = d.ci(d = d, n = df, alpha = .10)[1:2, 3],
        LL95 = d.ci(d = d, n = df)[1:2, 1],
        UL95 = d.ci(d = d, n = df)[1:2, 3]
    )
)

results3 <- bind_rows(results3, .id = "Measure")

results3 <- arrange(results3, contrast)

results3$Measure <- factor(results3$Measure, levels = levels(dta_diff$measure))

results3

A tibble: 18 × 13

Measure	contrast	n	B	SE	t	df	p	d	LL90	UL90	LL95	UL95
<fct>	<fct>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
FS	Handler v. Indirect and Direct	284	0.05268218	0.02880507	1.828920	280	6.847509e-02	0.2185978	0.02121749	0.41559079	-0.0165610806	0.45336320
PANAS-PA	Handler v. Indirect and Direct	284	0.04980548	0.04010202	1.241969	280	2.152872e-01	0.1484437	-0.04855939	0.34518044	-0.0862746357	0.38289513
SCS-R	Handler v. Indirect and Direct	284	0.08922721	0.03243980	2.750548	280	6.337236e-03	0.3287533	0.13053063	0.52638968	0.0926174239	0.56430906
SHS	Handler v. Indirect and Direct	284	0.10424736	0.03162582	3.296274	280	1.106110e-03	0.3939801	0.19512191	0.59214092	0.1570971344	0.63016870
Integration	Handler v. Indirect and Direct	283	0.21124955	0.05606994	3.767608	279	2.010893e-04	0.4511215	0.25126649	0.65017602	0.2130570060	0.68838876
Stress	Handler v. Indirect and Direct	284	-0.34808137	0.06948780	5.009244	280	9.690912e-07	-0.5987192	-0.79917741	-0.39721079	-0.8376882559	-0.35871912
Homesick	Handler v. Indirect and Direct	284	-0.10967163	0.06147159	1.784103	280	7.548976e-02	-0.2132411	-0.41020721	-0.01589082	-0.4479787796	0.02188227
Lonely	Handler v. Indirect and Direct	284	-0.06392091	0.01794412	3.562221	280	4.320507e-04	-0.4257668	-0.62421415	-0.22655941	-0.6623089552	-0.18847882
PANAS-NA	Handler v. Indirect and Direct	284	-0.13495201	0.03690386	3.656854	280	3.049329e-04	-0.4370776	-0.63563737	-0.23774431	-0.6737526733	-0.19963154
FS	Indirect v. Direct	284	0.02946717	0.02492737	1.182121	280	2.381606e-01	0.1412905	-0.05568076	0.33800799	-0.0933908061	0.37571909
PANAS-PA	Indirect v. Direct	284	0.07535830	0.03823624	1.970861	280	4.972458e-02	0.2355629	0.03806911	0.43263640	0.0002774757	0.47042926
SCS-R	Indirect v. Direct	284	0.06431532	0.03010435	2.136413	280	3.351274e-02	0.2553502	0.05772268	0.45252389	0.0199059268	0.49034192
SHS	Indirect v. Direct	284	0.06347497	0.02741215	2.315578	280	2.130405e-02	0.2767645	0.07897506	0.47405800	0.0411293311	0.51190496
Integration	Indirect v. Direct	283	0.12043551	0.04232461	2.845519	279	4.762586e-03	0.3407136	0.14202628	0.53879128	0.1040248954	0.57679777
Stress	Indirect v. Direct	284	-0.32432636	0.05512822	5.883128	280	1.146346e-08	-0.7031683	-0.90517639	-0.49993842	-0.9440059979	-0.46113593
Homesick	Indirect v. Direct	284	-0.15655633	0.05632442	2.779546	280	5.812134e-03	-0.3322193	-0.52988078	-0.13396324	-0.5678089787	-0.09604393
Lonely	Indirect v. Direct	284	-0.06668292	0.01587280	4.201080	280	3.574753e-05	-0.5021251	-0.70136930	-0.30199078	-0.7396352704	-0.26374064
PANAS-NA	Indirect v. Direct	284	-0.13813656	0.03263509	4.232762	280	3.132658e-05	-0.5059118	-0.70520037	-0.30572771	-0.7434746112	-0.26746881

Results of the ANCOVAs and Planned Contrasts

tbl3 <- results3

tbl3[, c("B", "SE", "t", "d", "LL95", "UL95")] <- sapply(
    tbl3[, c("B", "SE", "t", "d", "LL95", "UL95")],
    format_num
)

tbl3$p <- f_num(tbl3$p, digits = 3)
tbl3$p <- ifelse(tbl3$p == ".000", "<.001", tbl3$p)

tbl3$Measure <- as.character(tbl3$Measure)

keepers3 <- c("Measure", "B", "SE", "t", "p", "d", "LL95", "UL95")

tbl3 <- rbind(
    c("Handler v. Indirect and Direct", rep("", length(keepers3) - 1)),
    tbl3[tbl3$contrast == "Handler v. Indirect and Direct", keepers3],
    c("Indirect v. Direct", rep("", length(keepers3) - 1)),
    tbl3[tbl3$contrast == "Indirect v. Direct", keepers3]
)

write.csv(tbl3, row.names = FALSE, quote = FALSE)

Measure,B,SE,t,p,d,LL95,UL95
Handler v. Indirect and Direct,,,,,,,
FS,0.05,0.03,1.83,.068,0.22,−0.02,0.45
PANAS-PA,0.05,0.04,1.24,.215,0.15,−0.09,0.38
SCS-R,0.09,0.03,2.75,.006,0.33,0.09,0.56
SHS,0.10,0.03,3.30,.001,0.39,0.16,0.63
Integration,0.21,0.06,3.77,<.001,0.45,0.21,0.69
Stress,−0.35,0.07,5.01,<.001,−0.60,−0.84,−0.36
Homesick,−0.11,0.06,1.78,.075,−0.21,−0.45,0.02
Lonely,−0.06,0.02,3.56,<.001,−0.43,−0.66,−0.19
PANAS-NA,−0.13,0.04,3.66,<.001,−0.44,−0.67,−0.20
Indirect v. Direct,,,,,,,
FS,0.03,0.02,1.18,.238,0.14,−0.09,0.38
PANAS-PA,0.08,0.04,1.97,.050,0.24,0.00,0.47
SCS-R,0.06,0.03,2.14,.034,0.26,0.02,0.49
SHS,0.06,0.03,2.32,.021,0.28,0.04,0.51
Integration,0.12,0.04,2.85,.005,0.34,0.10,0.58
Stress,−0.32,0.06,5.88,<.001,−0.70,−0.94,−0.46
Homesick,−0.16,0.06,2.78,.006,−0.33,−0.57,−0.10
Lonely,−0.07,0.02,4.20,<.001,−0.50,−0.74,−0.26
PANAS-NA,−0.14,0.03,4.23,<.001,−0.51,−0.74,−0.27

ggplot(results3, aes(x = d, y = fct_rev(Measure))) +
    facet_wrap(~contrast, scales = "free_x") +
    geom_vline(xintercept = 0, size = 1.0, colour = "#e9e9e9") +
    geom_linerange(aes(xmin = LL95, xmax = UL95)) +
    geom_linerange(aes(xmin = LL90, xmax = UL90), size = 1) +
    geom_point() +
    theme_minimal(base_family = "Fira Sans") +
    labs(
        x = "Cohen's *d*",
        y = NULL
    ) +
    theme(
        axis.title.x = element_markdown(),
        strip.text = element_text(size = 12, hjust = 0, face = "bold")
    )

ggsave("h3h4.png", width = 6.5, height = 4, dpi = 320)

Exploratory Analyses

Controlling for Pet Ownership

Pet ownership was imbalanced between groups and could affect the outcome of the intervention. Are the results different if we control for pet ownership?

models_pet <- lapply(
    split(dta_diff, dta_diff$measure),
    function(x) aov(t2 ~ t1 + group + pets, data = x)
)

lapply(
    models_pet,
    function(x) lmtest::coeftest(x, vcov. = sandwich::vcovHC(x))
)

$FS

t test of coefficients:

                       Estimate Std. Error t value  Pr(>|t|)    
(Intercept)            0.645058   0.184039  3.5050 0.0005329 ***
t1                     0.896174   0.030691 29.2000 < 2.2e-16 ***
grouphandler_v_dog     0.050378   0.030058  1.6760 0.0948712 .  
groupindirect_v_direct 0.031426   0.025506  1.2321 0.2189636    
petsYes                0.030301   0.045047  0.6727 0.5017263    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


$`PANAS-PA`

t test of coefficients:

                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)             0.408241   0.170341  2.3966  0.01722 *  
t1                      0.883820   0.050764 17.4103  < 2e-16 ***
grouphandler_v_dog      0.053508   0.040843  1.3101  0.19125    
groupindirect_v_direct  0.067751   0.038819  1.7453  0.08205 .  
petsYes                -0.075160   0.071333 -1.0536  0.29297    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


$`SCS-R`

t test of coefficients:

                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            0.405041   0.221068  1.8322  0.06800 .  
t1                     0.931082   0.046986 19.8160  < 2e-16 ***
grouphandler_v_dog     0.082469   0.032976  2.5009  0.01297 *  
groupindirect_v_direct 0.063741   0.030959  2.0589  0.04044 *  
petsYes                0.056836   0.052674  1.0790  0.28153    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


$SHS

t test of coefficients:

                        Estimate Std. Error t value  Pr(>|t|)    
(Intercept)             0.791690   0.097189  8.1459 1.329e-14 ***
t1                      0.836609   0.026350 31.7502 < 2.2e-16 ***
grouphandler_v_dog      0.112726   0.031995  3.5233  0.000499 ***
groupindirect_v_direct  0.060847   0.027759  2.1920  0.029220 *  
petsYes                -0.091818   0.050611 -1.8142  0.070740 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


$Integration

t test of coefficients:

                        Estimate Std. Error t value  Pr(>|t|)    
(Intercept)             0.968152   0.132463  7.3088  2.96e-12 ***
t1                      0.750016   0.036877 20.3382 < 2.2e-16 ***
grouphandler_v_dog      0.229936   0.056431  4.0747  6.04e-05 ***
groupindirect_v_direct  0.121466   0.042863  2.8338  0.004941 ** 
petsYes                -0.114711   0.083601 -1.3721  0.171147    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


$Stress

t test of coefficients:

                        Estimate Std. Error t value  Pr(>|t|)    
(Intercept)             0.512033   0.156125  3.2796  0.001173 ** 
t1                      0.571372   0.049924 11.4449 < 2.2e-16 ***
grouphandler_v_dog     -0.354364   0.070238 -5.0452 8.244e-07 ***
groupindirect_v_direct -0.320181   0.055493 -5.7698 2.135e-08 ***
petsYes                 0.055079   0.105557  0.5218  0.602232    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


$Homesick

t test of coefficients:

                        Estimate Std. Error t value  Pr(>|t|)    
(Intercept)             0.478951   0.092890  5.1561 4.826e-07 ***
t1                      0.682987   0.043171 15.8204 < 2.2e-16 ***
grouphandler_v_dog     -0.130477   0.060764 -2.1473  0.032646 *  
groupindirect_v_direct -0.167069   0.055779 -2.9952  0.002993 ** 
petsYes                -0.047911   0.086584 -0.5533  0.580477    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


$Lonely

t test of coefficients:

                        Estimate Std. Error t value  Pr(>|t|)    
(Intercept)             0.050950   0.052964  0.9620  0.336906    
t1                      0.915644   0.025939 35.2992 < 2.2e-16 ***
grouphandler_v_dog     -0.060303   0.018180 -3.3170  0.001032 ** 
groupindirect_v_direct -0.066858   0.016286 -4.1052  5.33e-05 ***
petsYes                -0.012333   0.028176 -0.4377  0.661934    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


$`PANAS-NA`

t test of coefficients:

                        Estimate Std. Error t value  Pr(>|t|)    
(Intercept)             0.242868   0.101890  2.3836 0.0178221 *  
t1                      0.696194   0.055425 12.5610 < 2.2e-16 ***
grouphandler_v_dog     -0.136910   0.038350 -3.5700 0.0004212 ***
groupindirect_v_direct -0.135766   0.033189 -4.0907 5.653e-05 ***
petsYes                -0.023718   0.064705 -0.3666 0.7142325    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Composite Measure

Can we treat the ill-being and well-being items as measuring the same construct?

z <- dta

varnames <- c(
    "flourishing",
    "panas_pa",
    "scs",
    "shs",
    "integration",
    "stress",
    "homesickness",
    "loneliness",
    "panas_na"
)

keys <- c("stress", "homesickness", "loneliness", "panas_na")

pretty_varnames <- levels(dta_diff$measure)[1:9]

drop <- c("homesickness", "stress", "integration", "panas_pa")

ctt_t1 <- multi_ctt(
    x = z[z$time == "t1", ],
    varnames = varnames,
    keys = keys,
    drop = drop,
    pretty_varnames = levels(dta_diff$measure)[1:9]
)

write.csv(ctt_t1$rdrop, row.names = FALSE, quote = FALSE)
ctt_t1$alpha

Item,Step1,Step2,Step3,Step4,Step5
FS,.72,.73,.75,.80,.76
PANAS-PA,.42,.44,.44,.44,-
SCS-R,.78,.80,.83,.81,.83
SHS,.66,.68,.66,.70,.70
Integration,.34,.36,.39,-,-
Stress,.27,.28,-,-,-
Homesick,.22,-,-,-,-
Lonely,.75,.75,.77,.76,.78
PANAS-NA,.54,.56,.53,.55,.58

1. '.79'
2. '.82'
3. '.84'
4. '.87'
5. '.88'

ctt_t2 <- multi_ctt(
    x = z[z$time == "t2", ],
    varnames = varnames,
    keys = keys,
    drop = drop,
    pretty_varnames = levels(dta_diff$measure)[1:9]
)

write.csv(ctt_t2$rdrop, row.names = FALSE, quote = FALSE)
ctt_t2$alpha

Item,Step1,Step2,Step3,Step4,Step5
FS,.66,.69,.72,.76,.72
PANAS-PA,.34,.37,.38,.39,-
SCS-R,.75,.76,.81,.79,.80
SHS,.65,.66,.65,.68,.68
Integration,.39,.42,.43,-,-
Stress,.38,.36,-,-,-
Homesick,.23,-,-,-,-
Lonely,.78,.78,.80,.77,.80
PANAS-NA,.49,.49,.42,.42,.46

1. '.79'
2. '.82'
3. '.83'
4. '.84'
5. '.86'

Exploratory Factor Analysis

With all Variables

This is an EFA with all 9 variables, including the variables which the reliability analysis indicated were not well represented by the other variables (i.e., homesickness, stress, integration, and positive affect).

cormat <- cor(z[z$time == "t1", varnames], use = "pairwise.complete.obs")

Parallel Analysis

The number of factors according to parallel analysis is found by comparing the observed principal component eigenvalues to principal component eigenvalues of a sample random normal deviates with the same number of variables and observations.

The results below indicate that the first component explains more variance than random data eigenvalues, but the remaining components do not. As such, one factor should be retained.

set.seed(1234L)
cbind(
    ParallelAnalysis(n_obs = 284, n_vars = 9),
    observed = eigen(cormat, only.values = TRUE)$values
)

A matrix: 9 × 3 of type dbl

50%	95%	observed
1.2774524	1.3666655	4.10956694
1.1829021	1.2430673	1.04720639
1.1111807	1.1621184	0.96529336
1.0508972	1.0920646	0.85591454
0.9938312	1.0304687	0.71522799
0.9363919	0.9778409	0.50288798
0.8799228	0.9208879	0.44297754
0.8178225	0.8659311	0.28015716
0.7459072	0.8014340	0.08076809

Below shows the result of common factor analysis using maximum likelihood factor analysis.

ML1 is the loadings standardized loadings on the first (and only) factor (i.e., it is the correlation between the item and the factor). h2 is the communality. That is, the proportion of variance in the variable that is explained by the entire solution. Because it is a 1-factor solution, this is the same as the squared loading on the factor. u2 is just 1 - h2, or the unexplained variance.

Note that the soluation explains very little variance in homesickness, stress, integration, and positive affect.

fa_result <- fa(z[, varnames], nfactors = 1, fm = "ml")

fa_result

Factor Analysis using method =  ml
Call: fa(r = z[, varnames], nfactors = 1, fm = "ml")
Standardized loadings (pattern matrix) based upon correlation matrix
               ML1    h2    u2 com
flourishing   0.74 0.546 0.454   1
panas_pa      0.40 0.157 0.843   1
scs           0.95 0.909 0.091   1
shs           0.68 0.460 0.540   1
integration   0.48 0.232 0.768   1
stress       -0.30 0.089 0.911   1
homesickness -0.26 0.068 0.932   1
loneliness   -0.91 0.829 0.171   1
panas_na     -0.53 0.284 0.716   1

                ML1
SS loadings    3.57
Proportion Var 0.40

Mean item complexity =  1
Test of the hypothesis that 1 factor is sufficient.

The degrees of freedom for the null model are  36  and the objective function was  4.15 with Chi Square of  2339.5
The degrees of freedom for the model are 27  and the objective function was  0.59 

The root mean square of the residuals (RMSR) is  0.08 
The df corrected root mean square of the residuals is  0.09 

The harmonic number of observations is  568 with the empirical chi square  248.88  with prob <  9.1e-38 
The total number of observations was  568  with Likelihood Chi Square =  333.37  with prob <  1.5e-54 

Tucker Lewis Index of factoring reliability =  0.822
RMSEA index =  0.141  and the 90 % confidence intervals are  0.128 0.155
BIC =  162.13
Fit based upon off diagonal values = 0.96
Measures of factor score adequacy             
                                                   ML1
Correlation of (regression) scores with factors   0.97
Multiple R square of scores with factors          0.95
Minimum correlation of possible factor scores     0.89

With a Subset of Well-being Indicators

The reliability analysis and EFA both suggest that homesickness, stress, integration, and positive affect are not well represented by this common factor. Below shows an EFA with those items excluded.

Parallel Analysis

Parallel analysis shows even stronger evidence for a unidimensional solution. The observed eigenvalue for the second component is 0.58 (i.e., the second component explains 58% of the variance explained by one variable), which is much less than either the mean or 95th percentile random data eigenvalues.

set.seed(1234L)
cormat2 <- cor(
    z[z$time == "t1", varnames[-which(varnames %in% drop)]],
    use = "pairwise.complete.obs"
)

cbind(
    ParallelAnalysis(n_obs = 284, n_vars = 5),
    observed = eigen(cormat2, only.values = TRUE)$values
)

A matrix: 5 × 3 of type dbl

50%	95%	observed
1.1652737	1.2470997	3.47845623
1.0669686	1.1250723	0.58831769
0.9990498	1.0419693	0.52784112
0.9262531	0.9716274	0.32098205
0.8410315	0.9020703	0.08440292

EFA Results

The results of a maximum likelihood factor analysis with only flourishing, social connectedness, happiness, loneliness, and negative affect.

The loadings are all fairly large and the factor explains 60% of the variance in the 5 items.

fa_result2 <- fa(
    z[, varnames[-which(varnames %in% drop)]],
    nfactors = 1,
    fm = "ml"
)

fa_result2

Factor Analysis using method =  ml
Call: fa(r = z[, varnames[-which(varnames %in% drop)]], nfactors = 1, 
    fm = "ml")
Standardized loadings (pattern matrix) based upon correlation matrix
              ML1   h2    u2 com
flourishing  0.73 0.54 0.461   1
scs          0.96 0.93 0.074   1
shs          0.67 0.45 0.552   1
loneliness  -0.91 0.82 0.177   1
panas_na    -0.52 0.27 0.726   1

                ML1
SS loadings    3.01
Proportion Var 0.60

Mean item complexity =  1
Test of the hypothesis that 1 factor is sufficient.

The degrees of freedom for the null model are  10  and the objective function was  3.24 with Chi Square of  1827.29
The degrees of freedom for the model are 5  and the objective function was  0.24 

The root mean square of the residuals (RMSR) is  0.07 
The df corrected root mean square of the residuals is  0.1 

The harmonic number of observations is  568 with the empirical chi square  53.17  with prob <  3.1e-10 
The total number of observations was  568  with Likelihood Chi Square =  133.03  with prob <  5.4e-27 

Tucker Lewis Index of factoring reliability =  0.859
RMSEA index =  0.212  and the 90 % confidence intervals are  0.182 0.244
BIC =  101.31
Fit based upon off diagonal values = 0.99
Measures of factor score adequacy             
                                                   ML1
Correlation of (regression) scores with factors   0.98
Multiple R square of scores with factors          0.95
Minimum correlation of possible factor scores     0.90

Factor Scores

The factor scores are computed using regression with the standardized loadings.

Below is just wrangling the data to add each participants' factor score at each measurement occasion.

z$comp <- fa_result2$scores[, 1]

z <- z[, c("pid", "group", "time", "comp")]

z <- pivot_wider(
    z,
    names_from = "time",
    values_from = "comp"
)

z$measure <- "Composite"

z$diff21 <- z$t2 - z$t1

Hypothesis Tests with Factor Scores

We conducted the four null hypothesis significance tests using factor scores.

H1: Effect of the Intervention

Factor scores significantly increased pre-to-post intervention in all three conditions.

comp_h1 <- lapply(
    split(z, z$group),
    t_test, h = "post"
)

comp_h1 <- bind_rows(comp_h1, .id = "Condition")

rownames(comp_h1) <- NULL

comp_h1[, -1] <- round(comp_h1[, -1], 3)

comp_h1

A data.frame: 3 × 9

Condition	t	df	p	d	d90.LL	d90.UL	d95.LL	d95.UL
<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
Handler	1.862	93	0.033	0.192	0.020	0.363	-0.013	0.396
Indirect	3.796	94	0.000	0.389	0.213	0.564	0.180	0.597
Direct	7.767	94	0.000	0.797	0.601	0.989	0.564	1.026

H3 and H4: Effect of the Canines

Both contrasts were significant, indicating a significant effect of both the dog's presence, and direct physical contact with the canine.

contrasts(z$group) <- cbind(
    handler_v_dog = c(-1, 0.5, 0.5),
    indirect_v_direct = c(0, -1, 1)
)

model <- aov(t2 ~ t1 + group, data = z)
corrected_model <- lmtest::coeftest(model, vcov. = sandwich::vcovHC(model))

round(corrected_model, 3)

ds <- t2d(corrected_model[3:4, "t value"], n = 280)
round(d.ci(ds, n = c(280, 280)), 2)

t test of coefficients:

                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)               0.181      0.022   8.308   <2e-16 ***
t1                        0.935      0.039  24.253   <2e-16 ***
grouphandler_v_dog        0.116      0.031   3.778   <2e-16 ***
groupindirect_v_direct    0.093      0.030   3.133    0.002 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

A matrix: 2 × 3 of type dbl

	lower	effect	upper
grouphandler_v_dog	0.21	0.45	0.69
groupindirect_v_direct	0.14	0.37	0.61

dta$comp <- fa_result2$scores[, 1]

head(dta)

A data.frame: 6 × 14

	pid	group	time	pets	stress	homesickness	loneliness	panas_pa	panas_na	shs	scs	flourishing	integration	comp
	<int>	<fct>	<fct>	<chr>	<int>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<int>	<dbl>
1	1	Handler	t1	No	2	3	2.25	3.2	1.0	2.333333	3.900000	6.125	1	-0.5516554
2	1	Handler	t2	No	2	3	1.70	3.8	1.2	3.333333	3.800000	6.000	2	-0.3227439
3	2	Direct	t1	No	2	4	1.90	3.0	1.8	3.333333	5.150000	5.250	4	0.4879954
4	2	Direct	t2	No	1	2	1.60	3.2	1.0	4.000000	5.263158	6.000	4	0.8731664
5	3	Indirect	t1	No	4	3	2.25	2.8	1.6	2.666667	4.100000	5.375	3	-0.4837059
6	3	Indirect	t2	No	3	2	2.25	3.0	1.6	3.333333	4.150000	5.375	3	-0.4002878

cor_tbl_r <- lapply(split(dta[, c(varnames, "comp")], dta$time), function(x) round(corr.test(x)$r, 3))
cor_tbl_r

$t1

A matrix: 10 × 10 of type dbl

	flourishing	panas_pa	scs	shs	integration	stress	homesickness	loneliness	panas_na	comp
flourishing	1.000	0.524	0.732	0.684	0.282	-0.234	-0.213	-0.589	-0.498	0.763
panas_pa	0.524	1.000	0.394	0.363	0.230	-0.163	-0.074	-0.311	-0.197	0.405
scs	0.732	0.394	1.000	0.633	0.445	-0.230	-0.224	-0.889	-0.524	0.988
shs	0.684	0.363	0.633	1.000	0.263	-0.303	-0.196	-0.578	-0.468	0.690
integration	0.282	0.230	0.445	0.263	1.000	-0.028	-0.085	-0.419	-0.213	0.441
stress	-0.234	-0.163	-0.230	-0.303	-0.028	1.000	0.078	0.206	0.312	-0.250
homesickness	-0.213	-0.074	-0.224	-0.196	-0.085	0.078	1.000	0.277	0.125	-0.249
loneliness	-0.589	-0.311	-0.889	-0.578	-0.419	0.206	0.277	1.000	0.548	-0.931
panas_na	-0.498	-0.197	-0.524	-0.468	-0.213	0.312	0.125	0.548	1.000	-0.583
comp	0.763	0.405	0.988	0.690	0.441	-0.250	-0.249	-0.931	-0.583	1.000

$t2

A matrix: 10 × 10 of type dbl

	flourishing	panas_pa	scs	shs	integration	stress	homesickness	loneliness	panas_na	comp
flourishing	1.000	0.468	0.691	0.649	0.293	-0.244	-0.144	-0.625	-0.356	0.743
panas_pa	0.468	1.000	0.354	0.323	0.175	-0.163	-0.004	-0.308	-0.091	0.371
scs	0.691	0.354	1.000	0.611	0.480	-0.239	-0.212	-0.872	-0.432	0.986
shs	0.649	0.323	0.611	1.000	0.303	-0.311	-0.216	-0.591	-0.383	0.675
integration	0.293	0.175	0.480	0.303	1.000	-0.158	-0.064	-0.501	-0.216	0.492
stress	-0.244	-0.163	-0.239	-0.311	-0.158	1.000	0.198	0.295	0.445	-0.285
homesickness	-0.144	-0.004	-0.212	-0.216	-0.064	0.198	1.000	0.275	0.202	-0.241
loneliness	-0.625	-0.308	-0.872	-0.591	-0.501	0.295	0.275	1.000	0.444	-0.928
panas_na	-0.356	-0.091	-0.432	-0.383	-0.216	0.445	0.202	0.444	1.000	-0.482
comp	0.743	0.371	0.986	0.675	0.492	-0.285	-0.241	-0.928	-0.482	1.000

cor_tbl_p <- lapply(split(dta[, c(varnames, "comp")], dta$time), function(x) round(corr.test(x)$p, 3))
cor_tbl_p

$t1

A matrix: 10 × 10 of type dbl

	flourishing	panas_pa	scs	shs	integration	stress	homesickness	loneliness	panas_na	comp
flourishing	0	0.000	0	0.000	0.000	0.001	0.003	0.000	0.000	0
panas_pa	0	0.000	0	0.000	0.001	0.036	0.604	0.000	0.007	0
scs	0	0.000	0	0.000	0.000	0.001	0.002	0.000	0.000	0
shs	0	0.000	0	0.000	0.000	0.000	0.007	0.000	0.000	0
integration	0	0.000	0	0.000	0.000	0.637	0.604	0.000	0.003	0
stress	0	0.006	0	0.000	0.637	0.000	0.604	0.004	0.000	0
homesickness	0	0.211	0	0.001	0.151	0.190	0.000	0.000	0.175	0
loneliness	0	0.000	0	0.000	0.000	0.000	0.000	0.000	0.000	0
panas_na	0	0.001	0	0.000	0.000	0.000	0.035	0.000	0.000	0
comp	0	0.000	0	0.000	0.000	0.000	0.000	0.000	0.000	0

$t2

A matrix: 10 × 10 of type dbl

	flourishing	panas_pa	scs	shs	integration	stress	homesickness	loneliness	panas_na	comp
flourishing	0.000	0.000	0	0	0.000	0.000	0.061	0	0.000	0.000
panas_pa	0.000	0.000	0	0	0.022	0.035	0.952	0	0.379	0.000
scs	0.000	0.000	0	0	0.000	0.001	0.003	0	0.000	0.000
shs	0.000	0.000	0	0	0.000	0.000	0.003	0	0.000	0.000
integration	0.000	0.003	0	0	0.000	0.039	0.561	0	0.003	0.000
stress	0.000	0.006	0	0	0.008	0.000	0.006	0	0.000	0.000
homesickness	0.015	0.952	0	0	0.280	0.001	0.000	0	0.005	0.001
loneliness	0.000	0.000	0	0	0.000	0.000	0.000	0	0.000	0.000
panas_na	0.000	0.126	0	0	0.000	0.000	0.001	0	0.000	0.000
comp	0.000	0.000	0	0	0.000	0.000	0.000	0	0.000	0.000

cor_tbl_sig <- lapply(
    cor_tbl_p,
    function(x) {
        x <- ifelse(x < .001, "***", ifelse(x < .01, "**", ifelse(x < .05, "*", "")))
        x[upper.tri(x, diag = TRUE)] <- ""
        x
    }
)
cor_tbl_sig

$t1

A matrix: 10 × 10 of type chr

	flourishing	panas_pa	scs	shs	integration	stress	homesickness	loneliness	panas_na	comp
flourishing
panas_pa	***
scs	***	***
shs	***	***	***
integration	***	***	***	***
stress	***	**	***	***
homesickness	***		***	**
loneliness	***	***	***	***	***	***	***
panas_na	***	**	***	***	***	***	*	***
comp	***	***	***	***	***	***	***	***	***

$t2

A matrix: 10 × 10 of type chr

	flourishing	panas_pa	scs	shs	integration	stress	homesickness	loneliness	panas_na	comp
flourishing
panas_pa	***
scs	***	***
shs	***	***	***
integration	***	**	***	***
stress	***	**	***	***	**
homesickness	*		***	***		**
loneliness	***	***	***	***	***	***	***
panas_na	***		***	***	***	***	**	***
comp	***	***	***	***	***	***	***	***	***

cor_tbl_m <- lapply(split(dta[, c(varnames, "comp")], dta$time),
    colMeans,
    na.rm = TRUE
)
cor_tbl_m

cor_tbl_sd <- lapply(split(dta[, c(varnames, "comp")], dta$time),
    sapply,
    sd,
    na.rm = TRUE
)
cor_tbl_sd

$t1

flourishing

5.67787977867203

panas_pa

3.28908450704225

scs

4.41613232023721

shs

3.31103286384977

integration

3.22183098591549

stress

3.20774647887324

homesickness

2.08802816901408

loneliness

2.04018717568569

panas_na

1.97869718309859

comp

-0.0935518891776789

$t2

flourishing

5.74201458752515

panas_pa

3.2955985915493

scs

4.53423883946957

shs

3.53403755868545

integration

3.35335689045936

stress

2.37323943661972

homesickness

1.8943661971831

loneliness

1.91404744255004

panas_na

1.61285211267606

comp

0.0935518891776793

$t1

flourishing

0.78207518362171

panas_pa

0.682844192120855

scs

0.823891866633276

shs

0.875825679921248

integration

1.04834489760614

stress

0.970852522792361

homesickness

1.20481849984894

loneliness

0.530448811354172

panas_na

0.704838447163469

comp

0.960755754009622

$t2

flourishing

0.779902129442693

panas_pa

0.787192080670786

scs

0.866805525117641

shs

0.82273322631227

integration

1.02567104993462

stress

1.03034629966718

homesickness

1.1134116321464

loneliness

0.534739670867953

panas_na

0.671230055375456

comp

0.982357661533885

tbl5 <- lapply(cor_tbl_r, f_num, digits = 2)
tbl5 <- lapply(tbl5, function(x) gsub("-", "−", x))
tbl5 <- lapply(tbl5, matrix, ncol = 10)
tbl5 <- lapply(tbl5, function(x) {
    x[upper.tri(x, diag = TRUE)] <- ""
    x
})
tbl5$t1 <- matrix(paste0(tbl5$t1, cor_tbl_sig$t1), ncol = 10)
tbl5$t2 <- matrix(paste0(tbl5$t2, cor_tbl_sig$t2), ncol = 10)
tbl5 <- lapply(tbl5,
    matrix, ncol = 10,
    dimnames = list(
        c(levels(dta_diff$measure)[1:9], "Composite"),
        1:10
    )
)
tbl5$t1 <- cbind(
    Variable = rownames(tbl5$t1),
    M = format_num(cor_tbl_m$t1),
    SD = format_num(cor_tbl_sd$t1),
    tbl5$t1
)[, 1:12]
tbl5$t2 <- cbind(
    Variable = rownames(tbl5$t2),
    M = format_num(cor_tbl_m$t2),
    SD = format_num(cor_tbl_sd$t2),
    tbl5$t2
)[, 1:12]

tbl5 <- rbind(
    c("Pre-intervention", rep("", 11)),
    tbl5$t1,
    c("Post-intervention", rep("", 11)),
    tbl5$t2
)
rownames(tbl5) <- NULL
tbl5

A matrix: 22 × 12 of type chr

Variable	M	SD	1	2	3	4	5	6	7	8	9
Pre-intervention
FS	5.68	0.78
PANAS-PA	3.29	0.68	.52***
SCS-R	4.42	0.82	.73***	.39***
SHS	3.31	0.88	.68***	.36***	.63***
Integration	3.22	1.05	.28***	.23***	.44***	.26***
Stress	3.21	0.97	−.23***	−.16**	−.23***	−.30***	−.03
Homesick	2.09	1.20	−.21***	−.07	−.22***	−.20**	−.09	.08
Lonely	2.04	0.53	−.59***	−.31***	−.89***	−.58***	−.42***	.21***	.28***
PANAS-NA	1.98	0.70	−.50***	−.20**	−.52***	−.47***	−.21***	.31***	.12*	.55***
Composite	−0.09	0.96	.76***	.41***	.99***	.69***	.44***	−.25***	−.25***	−.93***	−.58***
Post-intervention
FS	5.74	0.78
PANAS-PA	3.30	0.79	.47***
SCS-R	4.53	0.87	.69***	.35***
SHS	3.53	0.82	.65***	.32***	.61***
Integration	3.35	1.03	.29***	.17**	.48***	.30***
Stress	2.37	1.03	−.24***	−.16**	−.24***	−.31***	−.16**
Homesick	1.89	1.11	−.14*	−.00	−.21***	−.22***	−.06	.20**
Lonely	1.91	0.53	−.62***	−.31***	−.87***	−.59***	−.50***	.30***	.28***
PANAS-NA	1.61	0.67	−.36***	−.09	−.43***	−.38***	−.22***	.44***	.20**	.44***
Composite	0.09	0.98	.74***	.37***	.99***	.68***	.49***	−.28***	−.24***	−.93***	−.48***