Replicating Table 1

---

This script replicates Table 1 of the manuscript (Descriptive statistics of Strava Clubs).

Getting started

clean up

rm (list = ls( ))

custom functions

• fpackage.check: Check if packages are installed (and install if not) in R (source)

fpackage.check <- function(packages) {
    lapply(packages, FUN = function(x) {
        if (!require(x, character.only = TRUE)) {
            install.packages(x, dependencies = TRUE)
            library(x, character.only = TRUE)
        }
    })
}

necessary packages

We install and load the packages we need later on: - RSiena - igraph: Descriptives (dyad/triad census, degrees) - dplyr: for data manipulation - sna: for network analysis - knitr: for generating tables - kableExtra: for manipulating tables

packages = c("RSiena", "igraph", "dplyr", "sna", "knitr", "kableExtra")
fpackage.check(packages)

additional functions

# density: observed relations divided by possible relations
fdensity <- function(x) {
    # x is your nomination network make sure diagonal cells are NA
    diag(x) <- NA
    # take care of RSiena structural zeros, set as missing.
    x[x == 10] <- NA
    sum(x == 1, na.rm = T)/(sum(x == 1 | x == 0, na.rm = T))
}

# calculate intragroup density
fdensityintra <- function(x, A) {
    # A is matrix indicating whether nodes in dyad have same node attributes
    diag(x) <- NA
    x[x == 10] <- NA
    diag(A) <- NA
    sum(x == 1 & A == 1, na.rm = T)/(sum((x == 1 | x == 0) & A == 1, na.rm = T))
}

# calculate intragroup density
fdensityinter <- function(x, A) {
    # A is matrix indicating whether nodes in dyad have same node attributes
    diag(x) <- NA
    x[x == 10] <- NA
    diag(A) <- NA
    sum(x == 1 & A != 1, na.rm = T)/(sum((x == 1 | x == 0) & A != 1, na.rm = T))
}

# construct dyad characteristic whether nodes are similar/homogenous
fhomomat <- function(x) {
    # x is a vector of node-covariate
    xmat <- matrix(x, nrow = length(x), ncol = length(x))
    xmatt <- t(xmat)
    xhomo <- xmat == xmatt
    return(xhomo)
}

# a function to calculate all valid dyads.
fndyads <- function(x) {
    diag(x) <- NA
    x[x == 10] <- NA
    (sum((x == 1 | x == 0), na.rm = T))
}

# a function to calculate all valid intragroupdyads.
fndyads2 <- function(x, A) {
    diag(x) <- NA
    x[x == 10] <- NA
    diag(A) <- NA
    (sum((x == 1 | x == 0) & A == 1, na.rm = T))
}


fscolnet <- function(network, ccovar) {
    # Calculate coleman on network level:
    # https://reader.elsevier.com/reader/sd/pii/S0378873314000239?token=A42F99FF6E2B750436DD2CB0DB7B1F41BDEC16052A45683C02644DAF88215A3379636B2AA197B65941D6373E9E2EE413
    
    fhomomat <- function(x) {
        xmat <- matrix(x, nrow = length(x), ncol = length(x))
        xmatt <- t(xmat)
        xhomo <- xmat == xmatt
        return(xhomo)
    }
    
    fsumintra <- function(x, A) {
        # A is matrix indicating whether nodes constituting dyad have same characteristics
        diag(x) <- NA
        x[x == 10] <- NA
        diag(A) <- NA
        sum(x == 1 & A == 1, na.rm = T)
    }
    
    # expecation w*=sum_g sum_i (ni((ng-1)/(N-1)))
    network[network == 10] <- NA
    ni <- rowSums(network, na.rm = T)
    ng <- NA
    for (i in 1:length(ccovar)) {
        ng[i] <- table(ccovar)[rownames(table(ccovar)) == ccovar[i]]
    }
    N <- length(ccovar)
    wexp <- sum(ni * ((ng - 1)/(N - 1)), na.rm = T)
    
    # wgg1 how many intragroup ties
    w <- fsumintra(network, fhomomat(ccovar))
    
    Scol_net <- ifelse(w >= wexp, (w - wexp)/(sum(ni, na.rm = T) - wexp), (w - wexp)/wexp)
    return(Scol_net)
}

fMoran.I <- function(x, weight, scaled = FALSE, na.rm = FALSE, alternative = "two.sided", rowstandardize = TRUE) {
    if (rowstandardize) {
        if (dim(weight)[1] != dim(weight)[2]) 
            stop("'weight' must be a square matrix")
        n <- length(x)
        if (dim(weight)[1] != n) 
            stop("'weight' must have as many rows as observations in 'x'")
        ei <- -1/(n - 1)
        nas <- is.na(x)
        if (any(nas)) {
            if (na.rm) {
                x <- x[!nas]
                n <- length(x)
                weight <- weight[!nas, !nas]
            } else {
                warning("'x' has missing values: maybe you wanted to set na.rm = TRUE?")
                return(list(observed = NA, expected = ei, sd = NA, p.value = NA))
            }
        }
        ROWSUM <- rowSums(weight)
        ROWSUM[ROWSUM == 0] <- 1
        weight <- weight/ROWSUM
        s <- sum(weight)
        m <- mean(x)
        y <- x - m
        cv <- sum(weight * y %o% y)
        v <- sum(y^2)
        obs <- (n/s) * (cv/v)
        if (scaled) {
            i.max <- (n/s) * (sd(rowSums(weight) * y)/sqrt(v/(n - 1)))
            obs <- obs/i.max
        }
        S1 <- 0.5 * sum((weight + t(weight))^2)
        S2 <- sum((apply(weight, 1, sum) + apply(weight, 2, sum))^2)
        s.sq <- s^2
        k <- (sum(y^4)/n)/(v/n)^2
        sdi <- sqrt((n * ((n^2 - 3 * n + 3) * S1 - n * S2 + 3 * s.sq) - k * (n * (n - 1) * S1 - 2 * n * 
            S2 + 6 * s.sq))/((n - 1) * (n - 2) * (n - 3) * s.sq) - 1/((n - 1)^2))
        alternative <- match.arg(alternative, c("two.sided", "less", "greater"))
        pv <- pnorm(obs, mean = ei, sd = sdi)
        if (alternative == "two.sided") 
            pv <- if (obs <= ei) 
                2 * pv else 2 * (1 - pv)
        if (alternative == "greater") 
            pv <- 1 - pv
        list(observed = obs, expected = ei, sd = sdi, p.value = pv)
    } else {
        if (dim(weight)[1] != dim(weight)[2]) 
            stop("'weight' must be a square matrix")
        n <- length(x)
        if (dim(weight)[1] != n) 
            stop("'weight' must have as many rows as observations in 'x'")
        ei <- -1/(n - 1)
        nas <- is.na(x)
        if (any(nas)) {
            if (na.rm) {
                x <- x[!nas]
                n <- length(x)
                weight <- weight[!nas, !nas]
            } else {
                warning("'x' has missing values: maybe you wanted to set na.rm = TRUE?")
                return(list(observed = NA, expected = ei, sd = NA, p.value = NA))
            }
        }
        # ROWSUM <- rowSums(weight) ROWSUM[ROWSUM == 0] <- 1 weight <- weight/ROWSUM
        s <- sum(weight)
        m <- mean(x)
        y <- x - m
        cv <- sum(weight * y %o% y)
        v <- sum(y^2)
        obs <- (n/s) * (cv/v)
        if (scaled) {
            i.max <- (n/s) * (sd(rowSums(weight) * y)/sqrt(v/(n - 1)))
            obs <- obs/i.max
        }
        S1 <- 0.5 * sum((weight + t(weight))^2)
        S2 <- sum((apply(weight, 1, sum) + apply(weight, 2, sum))^2)
        s.sq <- s^2
        k <- (sum(y^4)/n)/(v/n)^2
        sdi <- sqrt((n * ((n^2 - 3 * n + 3) * S1 - n * S2 + 3 * s.sq) - k * (n * (n - 1) * S1 - 2 * n * 
            S2 + 6 * s.sq))/((n - 1) * (n - 2) * (n - 3) * s.sq) - 1/((n - 1)^2))
        alternative <- match.arg(alternative, c("two.sided", "less", "greater"))
        pv <- pnorm(obs, mean = ei, sd = sdi)
        if (alternative == "two.sided") 
            pv <- if (obs <= ei) 
                2 * pv else 2 * (1 - pv)
        if (alternative == "greater") 
            pv <- 1 - pv
        list(observed = obs, expected = ei, sd = sdi, p.value = pv)
    }
    
    
}

load data objects

# load the RSiena objects
load("clubdata_rsiena_freq.RData")
load("clubdata_rsiena_vol.RData")
load("clubdata.RData")

---

Descriptives table

Replicate the descriptives table:

tab <- matrix(nrow=16, ncol=5)
rownames(tab) <- c("Density", "Average degree", "Reciprocity index", "Global transitivity index", "Jaccard index", "Gender segregation", "Network autocorrelation", "Direct ties", "Frequency (per week)", "Hours (per week)", "Distance decay", "Frequency (per week)", "Hours (per week)", "Weekly running frequency at t1", "Weekly running hours at t1", "Male sex (%)" )
colnames(tab) <- c("Club 1","Club 2", "Club 3", "Club 4", "Club 5")

for (i in 1:5) { # for each club
  # densities over time
  densities <- as.numeric(strsplit(substring(capture.output(clubdata_rsiena_freq[[i]])[12], 20), " ")[[1]])
  # mean (standard deviation)
  tab[1, i] <- paste(as.character(round(mean(densities), 3)), " (", as.character(round(sd(densities), 3)), ")", sep="")
}


# Average degrees
for (i in 1:5) { # for each club
  # calculate average out-degree over time
  d <- vector ()
  for (t in 1:12) {
    d[[t]] <- mean(igraph::degree(igraph::graph_from_adjacency_matrix(clubdata_rsiena_freq[[i]]$depvars$kudonet[,,t]), mode = "out"))
    
  }
  # mean (standard deviation)
  tab[2, i] <- paste(as.character(round(mean(d), 3)), " (", as.character(round(sd(d), 3)), ")", sep="")
}


# Reciprocity index
for (i in 1:5) { # for each club
  # calculate reciprocity index over time
  ri <- vector ()
  
  for (t in 1:12) {
    mut <- igraph::dyad.census(igraph::graph_from_adjacency_matrix(clubdata_rsiena_freq[[i]]$depvars$kudonet[,,t]))[[1]]
    nonmut <- igraph::dyad.census(igraph::graph_from_adjacency_matrix(clubdata_rsiena_freq[[i]]$depvars$kudonet[,,t]))[[2]]
    ri[[t]] <- mut/(mut+nonmut)
    
  }
  tab[3, i] <- paste(as.character(round(mean(ri), 3)), " (", as.character(round(sd(ri), 3)), ")", sep="")
}

# Transitivity index
for (i in 1:5) { # for each club
  # calculate transitivity index over time
  ti <- vector ()
  
  for (t in 1:12) {
    ti[[t]] <- igraph::transitivity(igraph::graph_from_adjacency_matrix(clubdata_rsiena_freq[[i]]$depvars$kudonet[,,t]))
    
  }
  tab[4, i] <- paste(as.character(round(mean(ti), 3)), " (", as.character(round(sd(ti), 3)), ")", sep="")
}

# Jaccard index

# Gender segregation 
for (i in 1:5) { # for each club
  # calculate coleman's homophily index over time
  col <- vector ()
  
  for (t in 1:12) {
    col[[t]] <- fscolnet(clubdata_rsiena_freq[[i]]$depvars$kudonet[,,t], clubdata_rsiena_freq[[i]]$cCovars$gender)
  }
  tab[6, i] <- paste(as.character(round(mean(col), 3)), " (", as.character(round(sd(col), 3)), ")", sep="")
}

# Network autocorrelation
# frequency
for (i in 1:5) { # for each club
  # we calculate moran's i
  # both with direct ties and indirect ties as well
  # for frequency and volume
  mor_freq_dir <- vector ()
  mor_freq_dist <- vector ()
  mor_vol_dir <- vector ()
  mor_vol_dist <- vector ()
  
  for (t in 1:12) {         
    freq <- clubdata[[i]]$freq_run[,,t]   # attributes at time t
    vol <- clubdata[[i]]$time_run[,,t]   
    knet <- network::as.network(clubdata_rsiena_freq[[i]]$depvars$kudonet[,,t]) # net a time t
    geodistances <- sna::geodist(knet, count.paths=T) # geodesic distance
    geodistances <- geodistances$gdist 
    diag(geodistances) <- Inf
    weights1 <- geodistances == 1               # only direct ties
    weights2 <- exp(-geodistances)              # also indirect
    
    mor_freq_dir[[t]] <- fMoran.I(freq, scaled = FALSE, # calculate moran on direct
                                  weight = weights1, 
                                  na.rm = TRUE, 
                                  rowstandardize = FALSE)[[1]] 
    mor_freq_dist[[t]] <- fMoran.I(freq, scaled = FALSE, # and indirect
                                   weight = weights2, 
                                   na.rm = TRUE, 
                                   rowstandardize = FALSE)[[1]]
    mor_vol_dir[[t]] <- fMoran.I(vol, scaled = FALSE, # calculate moran on direct
                                 weight = weights1, 
                                 na.rm = TRUE, 
                                 rowstandardize = FALSE)[[1]] 
    mor_vol_dist[[t]] <- fMoran.I(vol, scaled = FALSE, # and indirect
                                  weight = weights2, 
                                  na.rm = TRUE, 
                                  rowstandardize = FALSE)[[1]]
  }
  tab[9, i] <- paste(as.character(round(mean(mor_freq_dir), 3)), " (", as.character(round(sd(mor_freq_dir), 3)), ")", sep="")
  tab[10, i] <- paste(as.character(round(mean(mor_vol_dir), 3)), " (", as.character(round(sd(mor_vol_dir), 3)), ")", sep="")
  tab[12, i] <- paste(as.character(round(mean(mor_freq_dist), 3)), " (", as.character(round(sd(mor_freq_dist), 3)), ")", sep="")
  tab[13, i] <- paste(as.character(round(mean(mor_vol_dist), 3)), " (", as.character(round(sd(mor_vol_dist), 3)), ")", sep="")
}

# behavior at baseline (w1)
for (i in 1:5) { # for each club
  # running attributes at w1
  tab[14, i] <- paste(as.character(round(mean(clubdata[[i]]$freq_run[,,1], na.rm=T), 3)), " (", as.character(round(sd(clubdata[[i]]$freq_run[,,1], na.rm=T), 3)), ")", sep="")
  tab[15, i] <- paste(as.character(round(mean(clubdata[[i]]$time_run[,,1], na.rm=T), 3)), " (", as.character(round(sd(clubdata[[i]]$time_run[,,1], na.rm=T), 3)), ")", sep="")
}
# male sex (%)
for (i in 1:5) { # for each club
  # percentage male
  tab[16, i] <- round((length(clubdata[[i]][clubdata[[i]]$male==1])/clubdata[[i]]$netsize)*100, 3)
}

# for now we include jaccard manually; but we want to use a function to calculate these ideally.
tab[5, ] <- c("0.850 (0.060)", "0.75 (0.030)", "0.500 (0.060)", "0.800 (0.090)", "0.690 (0.020)")

Table

Table 1. Descriptive statistics of Strava clubs

	Club 1	Club 2	Club 3	Club 4	Club 5
Density	0.083 (0.014)	0.168 (0.014)	0.011 (0.002)	0.115 (0.018)	0.098 (0.013)
Average degree	2.392 (0.408)	10.156 (0.838)	1.787 (0.28)	1.385 (0.205)	7.418 (0.994)
Reciprocity index	0.803 (0.06)	0.634 (0.036)	0.386 (0.062)	0.72 (0.129)	0.607 (0.03)
Global transitivity index	0.581 (0.035)	0.544 (0.016)	0.276 (0.076)	0.735 (0.143)	0.526 (0.018)
Jaccard index	0.850 (0.060)	0.75 (0.030)	0.500 (0.060)	0.800 (0.090)	0.690 (0.020)
Gender segregation	0.235 (0.109)	0.294 (0.038)	0.093 (0.064)	-0.139 (0.12)	0.015 (0.034)
Network autocorrelation
Direct ties
Frequency (per week)	0.173 (0.169)	0.217 (0.081)	0.203 (0.069)	0.296 (0.264)	0.276 (0.04)
Hours (per week)	0.207 (0.151)	0.197 (0.099)	0.244 (0.065)	0.243 (0.302)	0.324 (0.054)
Distance decay
Frequency (per week)	0.108 (0.117)	0.112 (0.036)	0.141 (0.061)	0.311 (0.254)	0.157 (0.022)
Hours (per week)	0.138 (0.101)	0.107 (0.043)	0.161 (0.04)	0.275 (0.276)	0.18 (0.028)
Weekly running frequency at t1	2.033 (1.351)	3 (2.15)	1.885 (1.658)	1.385 (1.502)	3.208 (2.313)
Weekly running hours at t1	2.1 (1.561)	2.952 (2.161)	1.915 (1.726)	1.538 (1.613)	3.143 (2.475)
Male sex (%)	56.667	70.968	70.909	61.538	67.532
Note:
For network measures, mean values over time are given with standard deviations in parentheses. For behavioral attributes, means and standard deviations at baseline (t1: December) are given. Degrees reflect kudos ties, with awarding 1 or more kudos constituting the presence of a tie. The reciprocity index is the proportion of kudos ties that were reciprocated; the transitivity index is the number of closed triplets over the total number of triplets (both open and closed).
a The Jaccard similarity index measures the extent of tie change between consecutive waves
b The gender segregation measure used is Coleman’s homophily index.
c The network autocorrelation measure used is Moran’s I. We calculated Moran’s I for direct ties (i.e. with undirected path length 1) and for all ties to whom ego is (directly or indirectly) tied, using a distance-decay function for assigning weights. To construct the weight matrix, we measured the geodistance (d) for each dyad as the shortest (undirected) path length. We then used the negative exponential distance-decay function. We did not row-standardize.

---

Copyright © 2021 Rob Franken