nahodnaNP-noGraphs.R

## Toto je skript, který pomáhá určit, zda výsledky Q-metody mohly vzniknout náhodou, případně do jaké míry

## Vytvořil: 2020-10-22 FrK
## Upravil:  2020-11-05 FrK

## Encoding: windows-1250

## Hlavička
rm(list=ls())
dir = getwd()

library(dplyr)
library(qmethod)
library(foreign)
library(data.table)
library(stats)
library(here)
library(magrittr)
library(ggplot2)
library(ggpubr)



### Zcela náhodné vygenerování Q-sortů
## Definování potřebných funkcí 
# 23 stabilních
stb23 = function() {
  m = c(rep(1:48, 23)) %>% matrix(nrow = 48, byrow = FALSE) %>% data.frame()
  for (i in 1:23){
    m[,i] = sample(v)  
  }
  m
}

# 22 střídaných
str22 = function() {
  m = c(rep(1:48, 22)) %>% matrix(nrow = 48, byrow = FALSE) %>% data.frame()
  for (i in 1:22){
    m[,i] = sample(v)  
  }
  as.data.frame(m)
}

# Vzorek 28 resp. ze 45
vz28 = function(stab = b, stri = r) {
  df = cbind(stri[,sample(1:22, 5)], stab)
  names(df) = rNames
  df
}


## Definování potřebných konstant
# Pseudo-normální rozdělení
v = c(rep(-4, 3), rep(4, 3), rep(-3, 4), rep(3, 4), rep(-2, 6), rep(2, 6), rep(-1, 7), rep(1, 7), rep(0, 8))
length(v)  # test zda má 48 prvků
# Jména výroků
vNames = c(paste0("V0", 1:9), paste0("V", 10:48))

# Jména respondentů
rNames = c(paste0("R0", 1:9), paste0("R", 10:28))

# Parametry opakovaní
anals = 1000  # počet simulovaných analýz
samps = 1500  # počet vzorků pro jednu simulovanou analýzu

# Matice pro záznamy rozdílů
rozdilyPrumeruFaktoru  = c(rep(0, 48 * anals)) %>% matrix(ncol = 48)
prumerneRozdilyFaktoru = c(rep(0, 48 * anals)) %>% matrix(ncol = 48)


# Začátek simulace
start = Sys.time()

## Samotná simulace
for (analyza in 1:anals) {
  # Definice souboru
  b = stb23()
  r = str22()
  
  ## The first columns of two files
  # Loading data prepared in Stata
  results = qmethod(vz28(stab = b, stri = r), nfactors = 2, rotation = "cluster")[[6]]
  
  # Storing Factor1,2 from 1st file
  factor1 = results[[1]] %>% matrix (ncol = 48) %>% data.frame()
  factor2 = results[[2]] %>% matrix (ncol = 48) %>% data.frame()

  ## Loading all other columns
  if (samps>1) {
  
    for(n in 2:samps) {
      # Loading data
      results = qmethod(vz28(stab = b, stri = r), nfactors = 2, rotation = "cluster")[[6]]
   
      # Adding scores and factor1,2
      factor1 = rbind.data.frame(factor1, results[[1]])
      factor2 = rbind.data.frame(factor2, results[[2]])
    }
  }


  # Uložení výsledků do matice 'Průměrných rozdílů faktorů'
  names(factor1) = vNames
  names(factor2) = vNames
  rozdil = abs(factor1 - factor2) %>% gather("var", "v") %>% group_by(var) %>%  summarise(v = mean(v))
  prumerneRozdilyFaktoru[analyza, ] = rozdil$v %>% round(digits = 3)

  # Uložení výsledků do matice 'Rozdílů průměrnych faktorů'
  pf1 = factor1 %>% gather("var", "v") %>% group_by(var) %>%  summarise(v = mean(v)) %>% select(2)
  pf2 = factor2 %>% gather("var", "v") %>% group_by(var) %>%  summarise(v = mean(v)) %>% select(2)
  rozdilyPrumeruFaktoru[analyza, ] = abs(pf1$v - pf2$v) %>% round(digits = 3)

  # Počítadlo stavu
  for (variable in 1:500) {
    paste0("Toto byla analýza č. ", analyza, ".") %>% print()
  }
}

# Konec simulace
konec = Sys.time()
konec - start  # čas simulace

# Uložíme si výsledky simulace na později
saveRDS(rozdilyPrumeruFaktoru ,  file = 'RPF_2020-11-05.rds', ascii = TRUE)
saveRDS(prumerneRozdilyFaktoru , file = 'PRF_2020-11-05.rds', ascii = TRUE)


## Výpočty na základě rozdílů
# Příprava matice
dPRF = c(rep(0, 4 * anals)) %>% matrix(ncol = 4)
dRPF = c(rep(0, 4 * anals)) %>% matrix(ncol = 4)

# Spočteme a uložíme výskyty odlišností
for (i in 1:anals) {
  d[i,4] = (rozdily[i, ] > 4) %>% sum()  # Spočítáme, kolik "výroků" je odlišných extrémně,
  d[i,3] = (rozdily[i, ] <= 4 & rozdily[i, ] >= 3) %>% sum()  # kolik velmi,
  d[i,2] = (rozdily[i, ] < 3 & rozdily[i, ] > 1) %>% sum()  # kolik středně,
  d[i,1] = (rozdily[i, ] <= 1) %>% sum()  # kolik málo - všechny hranice jsou zvolené, aby odpovídali zpracování Q-grafu
}

# Nejprve kategorie jednotlivě
((d[,4] >= 05) %>% sum()) / anals * 100
((d[,4] <= 05) %>% sum()) / anals * 100
((d[,3] >= 12) %>% sum()) / anals * 100
((d[,3] <= 12) %>% sum()) / anals * 100
((d[,2] >= 10) %>% sum()) / anals * 100
((d[,2] <= 10) %>% sum()) / anals * 100
((d[,1] <= 21) %>% sum()) / anals * 100
((d[,1] >= 21) %>% sum()) / anals * 100
# A konečně zároveň!
(((d[,1] <= 21) & (d[,4] >= 5)) %>% sum()) / anals * 100
(((d[,1] >= 21) & (d[,4] >= 5)) %>% sum()) / anals * 100
(((d[,1] <= 21) & (d[,4] <= 5)) %>% sum()) / anals * 100
(((d[,1] >= 21) & (d[,4] <= 5)) %>% sum()) / anals * 100
(((d[,1] < 21) & (d[,4] > 5)) %>% sum()) / anals * 100
(((d[,1] > 21) & (d[,4] > 5)) %>% sum()) / anals * 100
(((d[,1] < 21) & (d[,4] < 5)) %>% sum()) / anals * 100
(((d[,1] > 21) & (d[,4] < 5)) %>% sum()) / anals * 100
# Výpis matice odlišností
d

# Histogramy s počty výroků v příslušných kategoriích odlišnosti
hist(d[,1])
hist(d[,2])
hist(d[,3])
hist(d[,4])

hist(rozdily)

# Kolik by rozdílů mělo být?
((rozdily > 4) %>% sum()) / anals
((rozdily <= 4 & rozdily >= 3) %>% sum()) / anals
((rozdily < 3 & rozdily > 1) %>% sum()) / anals
((rozdily <= 1) %>% sum()) / anals
mean(d[,1])
mean(d[,2])
mean(d[,3])
mean(d[,4])

# T-testem srovnáme, zda se simulace liší od našich výsledků
t.test(d[,1], alternative = "t", mu = 21)
t.test(d[,2], alternative = "t", mu = 10)
t.test(d[,3], alternative = "t", mu = 12)
t.test(d[,4], alternative = "t", mu = 5)





## Jak moc se budou lišit dva pohledy, když je vygenerujeme náhodně?
# Toto je distribuce pseudonormáního rozdělení pro Q-metodu
v = c(rep(-4, 3), rep(4, 3), rep(-3, 4), rep(3, 4), rep(-2, 6), rep(2, 6), rep(-1, 7), rep(1, 7), rep(0, 8))
length(v)  # test zda má 48 prvků

# Připravíme si matici pro numerickou simulaci - použijeme matici, protože je rychlejší ve výpočtech
rows = 1000000  # počet řádků matice
df = matrix(rep(0, rows * 4), ncol = 4)  # inicializace matice
colnames(df) = c('low', 'middle', 'high', 'vHigh')  # pojmenování sloupců matice

for (i in 1:rows) {
  x = sample(v)  # vytvoříme náhodnou kombinaci pseudonormálního rozdělení pro první pohled
  y = sample(v)  # totéž pro druhý
  d = abs(x - y)  # spočteme si rozdíly mezi pohledy v jednotlivých výrocích
  df[i,4] = (d > 4) %>% sum()  # Spočítáme, kolik "výroků" je odlišných extrémně,
  df[i,3] = (d <= 4 & d >= 3) %>% sum()  # kolik velmi,
  df[i,2] = (d < 3 & d > 1) %>% sum()  # kolik středně,
  df[i,1] = (d <= 1) %>% sum()  # kolik málo - všechny hranice jsou zvolené, aby odpovídali zpracování Q-grafu
}

# Histogramy s počty výroků v příslušných kategoriích odlišnosti
hist(df[,1])
hist(df[,2])
hist(df[,3])
hist(df[,4])

# Průměrné počty výroků v příslušných kategoriích odlišnosti
df %>% data.frame() %>% 
  summarise(low = mean(low), middle = mean(middle), high = mean(high), vHigh = mean(vHigh))

## Spočtení, kolik procent simulovaných distribucí vyšlo konsenzuálněji než ta naše
# Nejprve kategorie jednotlivě
((df[,4] <= 5) %>% sum()) / rows * 100
((df[,3] <= 12) %>% sum()) / rows * 100
((df[,2] >= 10) %>% sum()) / rows * 100
((df[,1] >= 21) %>% sum()) / rows * 100
# A konečně zároveň!
(((df[,1] >= 21) & (df[,4] <= 5)) %>% sum()) / rows * 100



## Další numerické simulace
# Jak je pravděpodobné, že z 24 náb. a 24 nenáb. výroků vybereme 5 náboženských?
n = c(rep(0, 24), rep(1, 24))  # Uděláme si vektor, 0 bude reprezentovat nenáb., 1 náb. výroky
length(n)  # ověření, že má vektor 48 prvků

# připravíme si vektor, kam budeme ukládat údaj, kolikrát jsme vybrali 5 náboženských výroků
pos = 1000000  # počet pozic, tj. délka vektoru
ex = rep(FALSE, pos)  # Připravíme prázdný vektor pro extrémní kontroverze
sl = rep(FALSE, pos)  # a pro shody dolní hranice CI
sh = rep(FALSE, pos)  # a pro shody, horní hranice CI
vl = rep(FALSE, pos)  # a pro vyhrocené shody, dolní
vh = rep(FALSE, pos)  # a pro vyhrocené shody, horní

# naplníme vektor
for (i in 1:pos) {
  ex[i] = (sample(n, 5) %>% sum())==5 
  sl[i] = (sample(n, 21) %>% sum())<8 
  sh[i] = (sample(n, 21) %>% sum())>13 
  vl[i] = (sample(n, 10) %>% sum())>7 
  vh[i] = (sample(n, 10) %>% sum())<3 
}

# Spočítáme procentní podíly a hranice konfidenčního intervalu
(ex %>% sum()) / pos * 100
(sl %>% sum()) / pos * 100
(sh %>% sum()) / pos * 100
(vl %>% sum()) / pos * 100
(vh %>% sum()) / pos * 100