# sets the directory of location of this script as the current directory
# setwd(dirname(rstudioapi::getSourceEditorContext()$path))
### load packages
require(pacman)
p_load('tidyverse', 'jsonlite', 'magrittr', 'xlsx',
'stargazer', 'psych', 'jtools', 'DT', 'ggstatsplot',
'lavaan',
'regsem', 'MplusAutomation', 'igraph', 'shiny', 'ggplot2', 'tidyLPA', 'MultilayerExtraction',
'Matrix', 'igraph', 'foreach', 'doParallel', 'parallel', 'R.matlab')
multilayer_extraction_results <- readRDS(file = "outputs/multilayer_extraction_results.Rda")
setwd("outputs/01_dataPreperation/final")
### load questionnaire
questionnaire <- readRDS(file = "questionnaire.rds")
CAMfiles <- readRDS(file = "CAMfiles.rds")
CAMdrawn <- readRDS(file = "CAMdrawn.rds")
CAMaggregated <- readRDS(file = "CAMaggregated.rds")
networkIndicators <- readRDS(file = "networkIndicators.rds")
CAMwordlist <- xlsx::read.xlsx2(file = "CAMwordlist.xlsx", sheetIndex = 1)
### load functions
# print(getwd())
setwd("../../../functions")
for(i in 1:length(dir())){
# print(dir()[i])
source(dir()[i], encoding = "utf-8")
}
setwd("../functions_CAMapp")
for(i in 1:length(dir())){
# print(dir()[i])
source(dir()[i], encoding = "utf-8")
}
rm(i)
### summary function
data_summary <- function(data, varname, groupnames){
require(plyr)
summary_func <- function(x, col){
c(mean = mean(x[[col]], na.rm=TRUE),
se = sd(x[[col]], na.rm=TRUE) / sqrt(length(x[[col]])))
}
data_sum<-ddply(data, groupnames, .fun=summary_func,
varname)
data_sum <- plyr::rename(data_sum, c("mean" = varname))
return(data_sum)
}Analyze CAM data on micro, meso, macro level for part II of basal attributes article
1 Notes
2 load data files
3 Micro level
Distributions, cluster analyses… on single node level…
3.1 Valence distributions
valenceData <- CAMfiles[[1]]
valenceData$value[valenceData$value == 10] <- 0
valenceDataShort <- valenceData %>% pivot_wider(names_from = text, values_from = value, id_cols = CAM)
hist1_names <- c("Akzeptanz eines neuen Materialsystems", "Energie generierend",
"Energie speichernd", "Insekten ähnlich", "aktive Formänderung durch Umwelteinwirkung",
"aktive Verhaltensänderung durch Umwelteinwirkung", "autonom",
"bioinspiriert", "elektronikfrei", "energieautonom", "energieeffizient",
"enthält Kunststoff", "haltbar", "intelligent", "langlebig",
"lebensähnlich")
hist2_names <- c("leicht zerstörbar", "multifunktional", "nachhaltig",
"passive Formänderung durch Umwelteinwirkung", "passive Verhaltensänderung durch Umwelteinwirkung",
"reaktionsfähig", "robust", "selbstheilend", "selbstreparierend",
"technologisch", "umweltfreundlich", "umweltschädlich", "wartungsfrei",
"wartungsintensiv", "widerstandsfähig", "zuverlässig", "ökologisch")
# adjust aesthetics of ggplot
ggplot_theme <- theme(axis.title.x = element_blank(),
axis.title.y = element_text(size=12),
axis.text.x = element_text(size=10,hjust=0.5,vjust=0.5,face="plain", colour = "black"),
axis.text.y = element_text(size=12,face="plain", colour = "black"),
panel.border = element_blank(),
axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
valenceData %>% filter(text %in% hist1_names) %>%
ggplot(aes(x = value)) + geom_histogram() + facet_wrap(~text, ncol = 4, nrow = 4) + ggplot_theme`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
valenceData %>% filter(text %in% hist2_names) %>%
ggplot(aes(x = value)) + geom_histogram() + facet_wrap(~text, ncol = 4, nrow = 5)`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# histograms sorted by MEAN VALENCE
hist1_names_mean <- CAMwordlist[order(CAMwordlist$mean_valence, decreasing = TRUE), "Words"][1:16]
hist2_names_mean <- CAMwordlist[order(CAMwordlist$mean_valence, decreasing = TRUE), "Words"][17:33]
valenceData %>% filter(text %in% hist1_names_mean) %>%
ggplot(aes(x = value)) + geom_histogram() + facet_wrap(~ factor(text, levels = hist1_names_mean), ncol = 4, nrow = 4) +
ggtitle("Sorted by avg. valence")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
valenceData %>% filter(text %in% hist2_names_mean) %>%
ggplot(aes(x = value)) + geom_histogram() + facet_wrap(~ factor(text, levels = hist2_names_mean), ncol = 4, nrow = 5) +
ggtitle("Sorted by avg. valence")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# histograms sorted by MEAN DEGREE
hist1_names_degree <- CAMwordlist[order(CAMwordlist$mean_degree, decreasing = TRUE), "Words"][1:16]
hist2_names_degree <- CAMwordlist[order(CAMwordlist$mean_degree, decreasing = TRUE), "Words"][17:33]
valenceData %>% filter(text %in% hist1_names_degree) %>%
ggplot(aes(x = value)) + geom_histogram() + facet_wrap(~ factor(text, levels = hist1_names_degree), ncol = 4, nrow = 4) +
ggtitle("Sorted by avg. degree")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
valenceData %>% filter(text %in% hist2_names_degree) %>%
ggplot(aes(x = value)) + geom_histogram() + facet_wrap(~ factor(text, levels = hist2_names_degree), ncol = 4, nrow = 5) +
ggtitle("Sorted by avg. degree")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# histograms sorted by MEAN TRANSITIVITY
hist1_names_transitivity <- CAMwordlist[order(CAMwordlist$mean_transitivity, decreasing = TRUE), "Words"][1:16]
hist2_names_transitivity <- CAMwordlist[order(CAMwordlist$mean_transitivity, decreasing = TRUE), "Words"][17:33]
valenceData %>% filter(text %in% hist1_names_transitivity) %>%
ggplot(aes(x = value)) + geom_histogram() + facet_wrap(~ factor(text, levels = hist1_names_transitivity), ncol = 4, nrow = 4) +
ggtitle("Sorted by avg. transitivity")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
valenceData %>% filter(text %in% hist2_names_transitivity) %>%
ggplot(aes(x = value)) + geom_histogram() + facet_wrap(~ factor(text, levels = hist2_names_transitivity), ncol = 4, nrow = 5) +
ggtitle("Sorted by avg. transitivity")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
4 Meso level
4.1 Multilayer community detection
4.1.1 MultilayerExtraction
Community detection with the MultilayerExtraction-Algorithm (https://github.com/jdwilson4/MultilayerExtraction)
setwd("outputs")
# recode vertex labels to integers
attributes_names <- c(hist1_names, hist2_names)
CAMaggregated[[4]][["node1"]] <- match(CAMaggregated[[4]][["node1"]], attributes_names)
CAMaggregated[[4]][["node2"]] <- match(CAMaggregated[[4]][["node2"]], attributes_names)
# recode layers (CAM IDs) to integers
layer_names <- unique(CAMfiles[[1]]$participantCAM)
CAMaggregated[[4]][["layer"]] <- match(CAMaggregated[[4]][["layer"]], layer_names)
# now we're ready for community extraction using the MultiLayer extraction package
set.seed(123)
start_time <- Sys.time()
#multilayer_extraction_results <- multilayer.extraction(adjacency = CAMaggregated[[4]], seed = 123, min.score = 0, prop.sample = .10)
end_time <- Sys.time()
end_time - start_timeTime difference of 0.0006108284 secs#saveRDS(multilayer_extraction_results, file = "multilayer_extraction_results.Rda")
plot(multilayer_extraction_results, main = "Diagnostic Plot")
# look at the first ... 20 communities
multilayer_extraction_communities <- data.frame(n_communities = seq(1:20))
multilayer_extraction_communities[paste0("community", 1:20)] <- NA
for(i in 1:20) {
object <- refine(multilayer_extraction_results, k = i, m = 192, n = nrow(CAMwordlist))
for(j in 1:i) {
multilayer_extraction_communities[i,j+1] <- paste(attributes_names[which(object$Vertices[,j] == 1)], collapse = ", ")
}
}
DT::datatable(multilayer_extraction_communities, options = list(pageLength = 5)) 4.1.2 GenLouvain
Community detection with the generalized Louvain algorithm as proposed by Mucha et al. (2010).
Re-import the results from GenLouvain as processed by MATLAB.
genlouvain <- readMat("outputs/genlouvain_results.mat", drop = "singletonLists")
genlouvain <- unlist(genlouvain, recursive = FALSE)
genlouvain_results <- list()
for(i in 1:length(genlouvain)) {
tmp_result <- unlist(unlist(genlouvain[[i]], recursive = FALSE)[[1]], recursive = FALSE)
tmp_result <- matrix(data = tmp_result, nrow = nrow(CAMwordlist), ncol = length(CAMaggregated[[5]]))
rownames(tmp_result) <- rownames(CAMaggregated[[5]][[1]])
colnames(tmp_result) <- names(CAMaggregated[[5]])
# named list with parameter 'omega' as name and attribute
attributes(tmp_result)$omega <- unlist(unlist(genlouvain[[i]], recursive = FALSE)[[4]], recursive = FALSE)
genlouvain_results[[paste0("omega", unlist(unlist(genlouvain[[i]], recursive = FALSE)[[4]], recursive = FALSE))]] <-
tmp_result
}Number of communities found:
for(i in 1:length(genlouvain_results)) {
print(paste("Found", length(unique(as.vector(genlouvain_results[[i]]))), "communities with ",
names(genlouvain_results)[i]))
}[1] "Found 953 communities with omega0"
[1] "Found 5 communities with omega0.1"
[1] "Found 5 communities with omega0.2"
[1] "Found 5 communities with omega0.3"
[1] "Found 4 communities with omega0.4"
[1] "Found 5 communities with omega0.5"
[1] "Found 5 communities with omega0.6"
[1] "Found 5 communities with omega0.7"
[1] "Found 5 communities with omega0.8"
[1] "Found 4 communities with omega0.9"
[1] "Found 5 communities with omega1"Take a look at the communities:
genlouvain_communities <- data.frame(omega = seq(0.1, 1, 0.1), community1 = NA, community2 = NA, community3 = NA,
community4 = NA, community5 = NA)
# omega = 0 with no coupling of the slices is not interesting, there are >900 communities.
# therefore we start at index 2 (omega = 0.1)
genlouvain_results <- genlouvain_results[-1]
for(i in 1:length(genlouvain_results)) {
for(j in 1:max(unlist(genlouvain_results))) {
genlouvain_communities[i, j+1] <- paste(rownames(genlouvain_results[[i]])[genlouvain_results[[i]][,1] == j], collapse=", ")
}
}
DT::datatable(genlouvain_communities, options = list(pageLength = 5)) 4.2 random graph - delete me
Compare the probability that two concepts are connected in randomly generated networks to drawn CAMs:
## get average number of drawn concepts (here fixed)
numConcepts <- mean(networkIndicators$num_nodes_macro)
## get average density
numDensity <- mean(networkIndicators$density_macro)
## simply get the average probability that two concepts are connected:
g <- igraph::random.graph.game(n = numConcepts, p.or.m = numDensity)
plot(g)
are.connected(g, 1, 2)Warning: `are.connected()` was deprecated in igraph 2.0.0.
ℹ Please use `are_adjacent()` instead.[1] FALSE#> whereby each edge to be drawn has the identical probability in the Erdős–Rényi model
vec_booleanConnected <- c()
# vec_booleanConnected2 <- c()
for(i in 1:10000){
g <- igraph::random.graph.game(n = numConcepts, p.or.m = numDensity)
vec_booleanConnected[i] <- are.connected(g, 1, 2)
# vec_booleanConnected2[i] <- are.connected(g, 22, 26)
}
baselineProbability <- mean(vec_booleanConnected)
baselineProbability[1] 0.0389# mean(vec_booleanConnected2)This baseline probability can be compared to all possible combinations of drawn concepts:
vec_boolean <- c()
for(i in 1:length(CAMdrawn)){
vec_boolean[i] <- are.connected(graph = CAMdrawn[[i]],
v1 = V(CAMdrawn[[i]])[V(CAMdrawn[[i]])$label == "robust"],
v2 = V(CAMdrawn[[i]])[V(CAMdrawn[[i]])$label == "widerstandsfähig"])
}
mean(vec_boolean)[1] 0.1720435 Macro level
5.1 plot aggregated CAM
# plot(CAMaggregated[[2]], vertex.size=diag(CAMaggregated[[1]]) / max(diag(CAMaggregated[[1]]))*20, edge.arrow.size=0.01)
# plot(CAMaggregated[[2]], vertex.size=(abs(V(CAMaggregated[[2]])$value)+1)*5, edge.arrow.size=0.01)
g = CAMaggregated[[2]]
g2 = simplify(CAMaggregated[[2]])
# plot(g2, edge.arrow.size=0.01,
# vertex.size=diag(CAMaggregated[[1]]) / max(diag(CAMaggregated[[1]]))*20)
E(g2)$weight = sapply(E(g2), function(e) {
length(all_shortest_paths(g, from=ends(g2, e)[1], to=ends(g2, e)[2])$res) } )
E(g2)$weight = E(g2)$weight * 2
# E(g2)$weight[E(g2)$weight == 1] <- NA
V(g2)$color[V(g2)$value <= .5 & V(g2)$value >= -.5] <- "yellow"
V(g2)$shape <- NA
V(g2)$shape <- ifelse(test = V(g2)$color == "yellow", yes = "square", no = "circle")
### > plot multiple times because of random layout
for(i in 1:3){
plot(g2, edge.arrow.size = 0,
layout=layout_nicely, vertex.frame.color="black", asp = .5, margin = -0.1,
vertex.size=diag(CAMaggregated[[1]]) / max(diag(CAMaggregated[[1]]))*5,
vertex.label.cex = .9,
edge.weight=2, edge.width=(E(g2)$weight/35))
}
CAMaggregated[[1]][1:5, 1:5] Akzeptanz eines neuen Materialsystems
Akzeptanz eines neuen Materialsystems 186
wartungsintensiv 49
bioinspiriert 51
intelligent 47
passive Formänderung durch Umwelteinwirkung 39
wartungsintensiv bioinspiriert
Akzeptanz eines neuen Materialsystems 31 24
wartungsintensiv 186 0
bioinspiriert 0 186
intelligent 2 6
passive Formänderung durch Umwelteinwirkung 3 3
intelligent
Akzeptanz eines neuen Materialsystems 39
wartungsintensiv 6
bioinspiriert 4
intelligent 186
passive Formänderung durch Umwelteinwirkung 10
passive Formänderung durch Umwelteinwirkung
Akzeptanz eines neuen Materialsystems 21
wartungsintensiv 1
bioinspiriert 4
intelligent 3
passive Formänderung durch Umwelteinwirkung 1866 Clusteranalyses
6.1 hierachical cluster analysis
hc_dat <- data.frame(CAM = unique(CAMfiles[[1]]$CAM), participantCAM = unique(CAMfiles[[1]]$participantCAM))
### create data set
## word vars
for(w in unique(CAMfiles[[1]]$text)){
varName_w <- str_remove_all(string = str_to_title(string = w, locale = "en"), pattern = " |\\W+")
hc_dat[[paste0("N_", varName_w)]] <- NA
hc_dat[[paste0("mean_", varName_w)]] <- NA
hc_dat[[paste0("SD_", varName_w)]] <- NA
}
## get N, mean, sd of single summarized concepts
verbose = FALSE
for(c in unique(CAMfiles[[1]]$CAM)){
if(verbose){
cat("considered CAM: ", c, "\n")
}
tmp_CAM_nodes <- CAMfiles[[1]][CAMfiles[[1]]$CAM == c, ]
tmp_CAM_nodes$value <- ifelse(test = tmp_CAM_nodes$value == 10, yes = 0, no = tmp_CAM_nodes$value)
for(w in unique(CAMfiles[[1]]$text)){
if(verbose){
cat("considered concept: ", w, "\n")
}
varName_w <- str_remove_all(string = str_to_title(string = w, locale = "en"), pattern = " |\\W+")
if(verbose){
cat(" > the freqeuncy, mean, SD are saved with the prefix N_, mean_, SD_ plus
word without white spaces: ", varName_w, "\n")
}
if(sum(tmp_CAM_nodes$text == w) > 0){
tmp_CAM_nodes_w <- tmp_CAM_nodes[tmp_CAM_nodes$text == w, ]
## add N
hc_dat[hc_dat$CAM == c, paste0("N_", varName_w)] <- nrow(tmp_CAM_nodes_w)
## add mean
hc_dat[hc_dat$CAM == c, paste0("mean_", varName_w)] <- mean(x = tmp_CAM_nodes_w$value)
## add SD, only if > 1
hc_dat[hc_dat$CAM == c, paste0("SD_", varName_w)] <- sd(x = tmp_CAM_nodes_w$value)
}
}
if(verbose){
cat("\n")
}
}
### compute clustering
hc_df <- hc_dat[, str_subset(string = colnames(hc_dat), pattern = "mean_")]
hc_df <- hc_df[, colSums(x = !is.na(hc_df)) >= 2] # only considers concepts which where drawn at least 2 times
hc_df_scaled <- scale(hc_df)
dist.eucl <- dist(hc_df_scaled, method = "euclidean")
hc_cluster <- hclust(dist.eucl, method = "ward.D2") # Ward's method
hc_cluster
Call:
hclust(d = dist.eucl, method = "ward.D2")
Cluster method : ward.D2
Distance : euclidean
Number of objects: 186 set_cutOff <- 20
plot(hc_cluster)
abline(h = set_cutOff, col = "red", lty = 2)
rect.hclust(hc_cluster, h=set_cutOff, border = "tomato")
groups <- cutree(hc_cluster, h=set_cutOff)
# groupsOut <- names(table(groups))[table(groups) >= 2]
table(groups)groups
1 2 3 4
53 47 72 14 aggregate(hc_df, by=list(cluster=groups), mean, na.rm = TRUE) cluster mean_AkzeptanzEinesNeuenMaterialsystems mean_Wartungsintensiv
1 1 0.39622642 -1.943396
2 2 0.21276596 -2.170213
3 3 0.45833333 -2.222222
4 4 0.07142857 -1.000000
mean_Bioinspiriert mean_Intelligent
1 1.03773585 1.3962264
2 1.29787234 2.0851064
3 0.43055556 0.8611111
4 0.07142857 1.0714286
mean_PassiveFormänderungDurchUmwelteinwirkung mean_Ökologisch
1 0.18867925 2.3396226
2 0.48936170 2.3404255
3 -0.50000000 1.6527778
4 -0.07142857 0.4285714
mean_LeichtZerstörbar mean_Technologisch mean_InsektenÄhnlich
1 -1.132075 0.9811321 0.1320755
2 -2.212766 0.7659574 0.1489362
3 -2.027778 0.7222222 -0.5000000
4 -1.357143 0.3571429 -0.2857143
mean_Wartungsfrei mean_Selbstheilend mean_Elektronikfrei
1 1.58490566 1.4716981 0.71698113
2 2.48936170 2.4255319 0.93617021
3 1.77777778 1.4166667 0.01388889
4 0.07142857 0.3571429 0.00000000
mean_Umweltfreundlich mean_Energieautonom mean_Selbstreparierend
1 2.849057 1.7358491 1.9056604
2 2.765957 1.7659574 2.5106383
3 2.402778 1.1527778 1.5000000
4 1.000000 0.3571429 0.2857143
mean_AktiveFormänderungDurchUmwelteinwirkung mean_Autonom
1 0.2830189 0.7169811
2 0.7872340 1.1063830
3 -0.2777778 0.5972222
4 -0.2857143 0.4285714
mean_Umweltschädlich mean_Multifunktional
1 -2.566038 1.8867925
2 -2.872340 2.0212766
3 -2.611111 1.4166667
4 -2.000000 0.9285714
mean_PassiveVerhaltensänderungDurchUmwelteinwirkung mean_Nachhaltig
1 0.18867925 2.377358
2 0.55319149 2.723404
3 -0.36111111 2.444444
4 0.07142857 1.000000
mean_AktiveVerhaltensänderungDurchUmwelteinwirkung mean_Robust
1 0.4716981 1.3207547
2 0.8510638 2.3404255
3 -0.1666667 2.0277778
4 -0.2857143 0.8571429
mean_Widerstandsfähig mean_Energieeffizient mean_Reaktionsfähig mean_Haltbar
1 1.2830189 2.1132075 0.7169811 1.43396226
2 2.4042553 2.4680851 1.1914894 2.38297872
3 1.8750000 2.1250000 0.3888889 2.29166667
4 0.4285714 0.7142857 0.2857143 0.07142857
mean_Zuverlässig mean_EnthältKunststoff mean_EnergieSpeichernd
1 2.0000000 -1.3396226 1.6415094
2 2.6170213 -1.1914894 1.9148936
3 2.3472222 -0.9444444 1.3750000
4 0.5714286 -0.5714286 0.4285714
mean_Lebensähnlich mean_Langlebig mean_EnergieGenerierend
1 0.52830189 1.4716981 1.7169811
2 0.74468085 2.8297872 2.2553191
3 -0.08333333 2.5416667 1.5694444
4 0.00000000 0.9285714 0.78571436.2 Latent profile analysis
just a preliminary analysis to get a feeling for what we’ve got.
# lpa <- valenceDataShort %>% select(-CAM) %>%
# estimate_profiles(1:9, variances = c("equal", "equal", "varying"),
# covariances = c("zero", "equal", "zero"))