WJS Ethics

Install required packages and set working directory

require(foreign)
require(FactoMineR)
require(sjmisc)
require(ade4)
require(factoextra)
require(soc.ca)
require(dplyr)
require(sjPlot)
require(car)
require(varhandle)
require(ggplot2)
require(plotluck)
require(beanplot)
require(seriation)
require(corrr)
require(paran)
require(tables)

Error in unloadNamespace(package) : 
  namnerommet <<'lattice'>> er importert av <<'coda', 'flexmix', 'Matrix', 'zoo', 'FactoMineR', 'lme4', 'psych', 'effects', 'sp', 'nlme'>>, og kan derfor ikkje lastast ut
Failed with error:  'Package 'lattice' version 0.20.34 cannot be unloaded'

require(xlsx)
require(corrplot)
require(ape)
require(grid)
require(matrixStats)
require(ggthemes)
require(data.table)
require(psych)
require(explor)
library(rworldmap)
library(CCA)

Import WJS consolidated data file and country data from SPSS

setwd("/Users/janfredrikhovden/Dropbox/DBXPAGAANDEARBEID/Statistikk/Rworkdir/WJS")
wjs<-read.spss("WJS Consolidated Data V4-01 310117.sav", to.data.frame=TRUE)
country<-read.spss("country data 100217.sav", to.data.frame=TRUE)
countryorg<-country
# wjs<-read.spss("WJS Consolidated Data V2-03 140916.sav", to.data.frame=TRUE)
# wjs<-read.spss("WJS Consolidated Data V3-01 261016.sav", to.data.frame=TRUE, use.missings=FALSE) # If one want to include all missing values
## NOTE: IF THE NUMBER OR ORDER OF COUNTRIES ARE CHANGED, CHECK COUNTRY ROW NUMBER belows where #CHECKROWNUMBER are added
wjs$ID<-as.numeric(rownames(wjs)) # assign row number as ID
wjc<-dplyr::left_join(wjs,country, by="COUNTRY")

Filter additional countries and make subsamples if needed

Functions

scale_by_rank <- function(x) apply(x, 2, rank)
Mean <- function(x) base::mean(x, na.rm=TRUE) # Mean ignoring missing values
Sd<-function(x) base::sd(x, na.rm=TRUE) # SD ignoring missing values
delete.na <- function(DF, n=0) {
  DF[rowSums(is.na(DF)) <= n,]
} # drop rows with a treshold for NAs

Recoding

# fixing naming of United Arab Emirates (should NOT be United Arabic Emirates)
levels(country$COUNTRY)[match("UAE",levels(country$COUNTRY))] <- "United Arab Emirates"
levels(wjs$COUNTRY)[match("UAE",levels(wjs$COUNTRY))] <- "United Arab Emirates"
# import and selection of variables for analysis
indep<- dplyr::select(wjs, COUNTRY, C1, C2, C6, C7, C20, C14A:C14K, C22:C23, T5, T7.1:T9, AUTONOMY, ID)
#recoding where necessary
#C1 manager
indep$manager<-recode(indep$C1, "c('Reporter', 'News writer', 'Trainee', 'Other')='not.manager'; NA=NA; else='manager'")
indep$manager <- factor(indep$manager, levels=c("not.manager", "manager")) # changer order of levels
indep$T5<-factor(indep$T5, levels=c("Rank-and-file","\"Junior\" manager","Senior/executive manager"))
#C2 fulltime
indep$fulltime<-recode(indep$C2, "c('Full-time employment')='Fulltime'; NA=NA; else='not.Fulltime'")
#C6/7 newsbeat
indep$hardnews<-recode(indep$C7, "c('News/current affairs', 'Politics', 'Foreign politics', 'Domestic politics', 'Economy', 'Crime & law', NA)='hardnews or general'; else='softnews'")
#C20 education
indep$highedu<-recode(indep$C20, "c('Not completed high school', 'Completed high school', 'some university studies, but no degree')='no.highedu'; NA=NA; else='high.edu'")
indep$master<-recode(indep$C20, "c('Not completed high school', 'Completed high school', 'some university studies, but no degree', 'College/Bachelor\342\200\231s degree or equivalent')='no.master';NA=NA; else='master'")
indep$edu3<-recode(indep$C20, "c('Not completed high school', 'Completed high school')='edu.nohigh';c('some university studies, but no degree', 'College/Bachelor\342\200\231s degree or equivalent')='edu.mid'; NA=NA; else='edu.master'")
indep$edu3 <- factor(indep$edu3, levels=c("edu.nohigh", "edu.mid", "edu.master")) # changer order of levels
# T7 medium (how?)
indep$broadcast<-recode(indep$T7.4, "c('Yes')='Broadcasting'; NA=NA; else='not.broadcasting'")
indep$broadcast[indep$T7.5=="Yes"]<-"Broadcasting"
indep$broadcast<-factor(indep$broadcast, levels=c("not.broadcasting", "Broadcasting"))
# T8 national media
indep$national<-recode(indep$T8, "c('National')='National'; NA=NA; else='not.national'")
indep$national<-factor(indep$national, levels=c("not.national", "National"))
# T9 ownership
indep$owner<-recode(indep$T9, "c('Purely public ownership', 'Purely state ownership', 'Mixed ownership but mostly public', 'Mixed ownership but mostly state-owned')='public/state'; NA=NA; else='private'")

C13 ETHICS

## Subsetting dataset and renaming variables
BAT<-dplyr::select(wjs,COUNTRY, ID, C13A:C14K)
# BAT<-BAT %>% filter(complete.cases(.))  # note filter on only complete cases
BAT<-plyr::rename(BAT, c(C13A="Always.follow.codes", C13B="Depends.on.situation", C13C="Matter.of.personal.judgement", C13D="Can.be.set.aside",C14A="Paying.for.info", C14B="Using.conf.documents", C14C="False.identity", C14D="Pressure.informants",
C14E="Use.personal.doc", C14F="Undercover.empl", C14G="Hidden.miccam", C14H="Recreations", C14J="Publ.unverified", C14K="Accept.money"))

#C13
BATm <- apply(dplyr::select(BAT, Always.follow.codes:Can.be.set.aside), 2, function(x) {x <- recode(x,"'strongly disagree'=1; 'somewhat disagree'=2; 'undecided'=3; 'somewhat agree'=4; 'strongly agree'=5"); x})
### re-add country and add labels
BATm<-data.frame(BATm)
BATm$COUNTRY<-BAT$COUNTRY 
BATm$ID<-BAT$ID 
BATmC13<-BATm # make a copy of the dataset (orginal)
BATmCOMPLETE<-BATm %>% filter(complete.cases(.)) # note filter on only complete C13* cases
BATmORG<-BATm # make a copy of the dataset (orginal)
### convert to matrix of averages by country and add rownames and labels
BATm<-BATm %>%
  group_by(COUNTRY) %>%
  summarise_each(funs(Mean), Always.follow.codes:Can.be.set.aside)
BATm<-data.frame(BATm)
rownames(BATm)<-BATm[,1]
BATmORGC13map<-BATm
BATm<-BATm[,2:5]
BATmORGC13<-BATm
#BATm<-BATm[c(-61),]# remove Tanzania 
### make table of mean and sd by question and country
DT <- data.table(BATmC13)
wide <- setnames(DT[, sapply(.SD, function(x) list(n=sum(complete.cases(x)), mean=round(mean(x, na.rm=TRUE), 2), sd=round(sd(x, na.rm=TRUE), 2))), by=COUNTRY], c("COUNTRY", sapply(names(DT)[-6], paste0, c(".n", ".men", ".SD"))))
table.meansdC13<-as.data.frame(wide[,1:13])

Maps C13

sPDF<-joinCountryData2Map(BATmORGC13map, joinCode = "NAME", nameJoinColumn = "COUNTRY")

67 codes from your data successfully matched countries in the map
0 codes from your data failed to match with a country code in the map
176 codes from the map weren't represented in your data

par(mai=c(0,0,0.2,0),xaxs="i",yaxs="i")
mapCountryData(sPDF, nameColumnToPlot = "Always.follow.codes", colourPalette  =c('black','grey90'))

2 colours specified and 7 required, using interpolation to calculate colours

mapCountryData(sPDF, nameColumnToPlot = "Depends.on.situation", colourPalette  =c('black','grey90'))

2 colours specified and 7 required, using interpolation to calculate colours

mapCountryData(sPDF, nameColumnToPlot = "Matter.of.personal.judgement", colourPalette  =c('black','grey90'))

2 colours specified and 7 required, using interpolation to calculate colours

mapCountryData(sPDF, nameColumnToPlot = "Can.be.set.aside", colourPalette  =c('black','grey90'))

2 colours specified and 7 required, using interpolation to calculate colours

Various descriptive statistics and figures C13

# sjp.setTheme(geom.outline.size = 0, geom.label.size = 3, title.size = 1.5)
#sjp.likert(BAT[,2:5], sort.frq = "pos.asc", show.prc.sign = TRUE, digits=0, show.n=FALSE, show.legend=FALSE, catcount=5,title="Agreement?     (C13, all countries)") # lickert plot
# Average mean and spread of roles (based on country averages), boxplot and beanplot
sortedBATm<-BATm[ , order(colMeans(BATm))] # sort by column mean
par(mar = c(4,15,2,2)) # set margins
boxplot(x = as.list(sortedBATm), horizontal = TRUE, las=1, main="Averages aggregated by country (C13)", log="")
par(mar = c(4,15,2,2)) # set margins

beanplot(x = as.list(sortedBATm), horizontal = TRUE, las=1, main="Averages aggregated by country (C13)", log="")

##Bertin plot (centered)
cBATm<-t(BATm)
cBATm<-scale(cBATm, scale=FALSE) #centering of means
cBATm<-t(cBATm)
cBATmC13<-cBATm
x <- scale_by_rank(as.matrix(scale(cBATm, scale=FALSE)))
order <- seriate(x, method="PCA")
bertinplot(x, order, options = list(panel=panel.squares, shading=TRUE, reverse=TRUE, mar = c(1,1,10,6), gp_labels=gpar(fontsize=8)))

order <- seriate(x, method="PCA_angle")
bertinplot(x, order, options = list(panel=panel.squares, shading=TRUE, reverse=TRUE, mar = c(1,1,10,6), gp_labels=gpar(fontsize=8)))

heatmap(cBATm, margins=c(5,5)) # heatmap (centered)

# network plot of correlations
c<-correlate(BATm)
network_plot(c, min_cor=.2, colors=c("red", "green")) # not working for now

c<-cor(BATm)
corrplot(c, type = "upper", order = "FPC", tl.col = "black", tl.srt = 45, diag=FALSE, tl.cex=0.8,mar=c(1,0,0,1))

#treemap
#library(treemap)
#tm <- treemap(
#  wjs,
#  index=c("C13A","C22"),
#  vSize="C23",
#  vColor="C23",
#  type="value"
#  )
# some nice HTML tables
#sjt.frq(BATm)
#sjt.corr(BATm, triangle="lower", digits=2)

PCA and biplot, mean centered

#BATm2<- BATm[-c(22,54), ] #Removal Ethiopia and Singapore as outliers
BATm2<- BATm
cBATm<-t(BATm2)
cBATm<-scale(cBATm, scale=FALSE) #centering of means
cBATm<-t(cBATm)
x<-t(cBATm)
x<-scale(x, scale=FALSE)
x<-t(x)
C13x<-x
paran(x)

Using eigendecomposition of correlation matrix.
Computing: 10%  20%  30%  40%  50%  60%  70%  80%  90%  100%


Results of Horn's Parallel Analysis for component retention
120 iterations, using the mean estimate

-------------------------------------------------- 
Component   Adjusted    Unadjusted    Estimated 
            Eigenvalue  Eigenvalue    Bias 
-------------------------------------------------- 
1           1.578132    1.866935      0.288803
2           1.112951    1.191158      0.078207
3           1.034434    0.941905     -0.09252
-------------------------------------------------- 

Adjusted eigenvalues > 1 indicate dimensions to retain.
(3 components retained)

res.pca<-PCA(x , scale.unit=FALSE, ncp=3, graph = TRUE)

#res.pca$var$coord
#sjt.pca(x)
fviz_pca_biplot(res.pca, geom = "text",  title="Ethics(C13)", axes =c(1,2), labelsize=2.5, repel=TRUE) + theme_minimal()

fviz_screeplot(res.pca, ncp=10)

fviz_contrib(res.pca, choice = "var", axes = 1)

fviz_contrib(res.pca, choice = "var", axes = 2)

fviz_contrib(res.pca, choice = "var", axes = 3)

indroles <- get_pca_ind(res.pca)
C13ind<-indroles$coord
# HCPC Clustering
res.pca.hcpc<-HCPC(res.pca, nb.clust=-1) # 4 clusters

# fviz_cluster(res.pca.hcpc, data = x, ellipse.type = "convex", repel=TRUE)+ theme_tufte() + labs(title = "C13 clusters (after PCA)") + ylab(" <--- 2: Absolutism") + xlab("1: <--- Subjectivism") + theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold"))
#res.pca.hcpc$desc.var
#K-means
fviz_nbclust(x, kmeans, method = "gap_stat") + labs(title = "C13 clusters (K-means)") # 1 cluster 

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 100)  [one "." per sample]:
.................................................. 50 
.................................................. 100 

# Non-metric MDS
library(ggrepel)
d<-dist(C13x[,c(-2)]) # also removing Argentina #CHECKROWNUMBER
fit <- isoMDS(d, k=2)

initial  value 6.020657 
iter   5 value 4.966264
final  value 4.895739 
converged

fit

$points
                             [,1]         [,2]
Albania               0.258984409  0.011497744
Argentina            -0.183570909 -0.210079535
Australia            -0.393562214  0.176350278
Austria              -0.623322790 -0.202296997
Bangladesh           -0.205624441 -0.432094209
Belgium              -0.008643381  0.236230278
Bhutan                0.402496417 -0.031589803
Botswana             -0.362219380  0.363285978
Brazil               -0.269438860  0.051660734
Bulgaria             -0.360805245  0.424142264
Canada               -0.175686078  0.351887873
Chile                -0.260635740  0.026905835
China                -0.216437701  0.149165365
Colombia             -0.354620923 -0.523853565
Croatia              -0.505291259 -0.362948555
Cyprus               -0.251302049  0.429038206
Czech Republic        0.500637945  0.534631530
Denmark               0.541397792  0.348185950
Ecuador              -0.207314845 -0.353619630
Egypt                 0.080423057 -0.570580364
El Salvador          -0.220828297 -0.198866658
Estonia              -0.319421736  0.032870988
Ethiopia              1.518482575  0.450973682
Finland              -0.583692325 -0.285093926
France                0.145684302  0.295620757
Germany              -0.640942868 -0.019179919
Greece               -0.141870588  0.279220327
Hong Kong             0.313294915 -0.174626607
Hungary               0.195409661  0.224304136
Iceland              -0.093826881  0.003848489
India                 0.610349564 -0.035917518
Indonesia            -0.504426247  0.178738686
Ireland               0.123165523  0.044182000
Israel                0.066644785  0.243173097
Italy                -0.571513391  0.321816975
Japan                 0.385008044  0.043034576
Kenya                 0.059837813  0.144617944
Kosovo               -0.280403006  0.077302431
Latvia               -0.060641746 -0.158310413
Malawi               -0.090198382  0.288509214
Malaysia              0.430909299 -0.212007973
Mexico                0.003937809 -0.095850547
Moldova              -0.150319640  0.085610526
Netherlands           0.469611444  0.005019157
New Zealand          -0.177517115 -0.020797823
Norway               -0.181271007 -0.233866555
Oman                  0.837770920 -0.488560733
Philippines           0.042946152 -0.194060550
Portugal             -0.484417685  0.295831773
Qatar                 0.386792796 -0.598292965
Romania              -0.294352932  0.076140633
Russia                0.467323649 -0.059511023
Serbia               -0.614377117 -0.321198665
Sierra Leone          0.330470860  0.441689895
Singapore             1.344721270  0.032266233
South Africa         -0.251214011  0.074651403
South Korea          -0.300339219 -0.146548562
Spain                -0.363942061  0.009376137
Sudan                 0.781981584 -0.637407781
Sweden                0.397166585  0.345554651
Switzerland          -0.344285506  0.131001730
Tanzania              0.641630160 -0.144360929
Thailand              0.739946406 -0.155975495
Turkey               -0.313832796 -0.180811402
United Arab Emirates -0.069613572 -0.336509578
UK                    0.020145858 -0.004417817
USA                  -0.665447649  0.160898620

$stress
[1] 4.895739

fit.sh<-Shepard(d, fit$points)
plot(fit.sh, pch = ".")
lines(fit.sh$x, fit.sh$yf, type = "S")

x <- fit$points[,1]
y <- fit$points[,2]
plot(x, y, xlab="Coordinate 1", ylab="Coordinate 2",
main="NonMetric MDS",   type="n")
text(x, y, labels = row.names(BATm), cex=.7)

space<-as.data.frame(fit$points)
space$label<-C13x[,-1]
qplot(V1, V2, data=space) + geom_text_repel(aes(label=row.names(fit$points))) + ggtitle("Non-parametric MDS C13") + theme_minimal()

#res.mds.hcpc<-HCPC(space, nb.clust=-1) 
#fviz_cluster(res.mds.hcpc, data = x, frame.type = "convex", repel=TRUE, labelsize=3)+ theme_tufte() + labs(title = "C13 clusters (after MDS)") 
# res.mds.hcpc$desc.var

PCA ETHICS C13 (non-aggregated)

ethic<-BATmC13
#roles<-dplyr::filter(roles,COUNTRY!='Singapore_NON-REP' & COUNTRY !='Ethiopia') # remove Singapore and Ethiopia
ethic<-ethic %>% filter(complete.cases(.))

x<-t(ethic[,c(-(5:6))])
x<-scale(x, scale=FALSE)
x<-t(x)
center_ethic<-as.data.frame(x)
center_ethic$COUNTRY<-ethic$COUNTRY
center_ethic$ID<-ethic$ID

KMO(as.matrix(center_ethic[,c(-(5:6))]))
paran(as.matrix(center_ethic[,c-(5:6)])) #only 1 component, but 3 with centered roles

# non-rotated solution, with centered data
res.pca<-PCA(center_ethic, quali.sup=5, quanti.sup = 6, scale.unit=TRUE, ncp=3, graph = FALSE) 
res.pca$var
summary.PCA(res.pca)
sweep(res.pca$var$coord,2,sqrt(res.pca$eig[1:ncol(res.pca$var$coord),1]),FUN="/")
#save coordinates
center_ethic$dim<-res.pca$ind$coord

fviz_screeplot(res.pca, ncp=10)
fviz_contrib(res.pca, choice = "var", axes = 1)
fviz_contrib(res.pca, choice = "var", axes = 2)
fviz_contrib(res.pca, choice = "var", axes = 3)

fviz_pca_biplot(res.pca, geom = "text",  invisible="ind", title="Ethics(C13)", axes =c(1,2), repel=TRUE) + theme_minimal()
fviz_pca_biplot(res.pca, geom = "text",  invisible="ind", title="Ethics(C13)", axes =c(1,3), repel=TRUE) + theme_minimal()

quali.coord <- as.data.frame(res.pca$quali.sup$coord)
quali.coord$name <- rownames(quali.coord)
quali.coord<-transform(quali.coord, Dim.1=Dim.1, Dim.2=Dim.2, Dim.3=Dim.3)
fviz_pca_var(res.pca, 
             invisible=c("ind","var", "quanti.sup"), axes=c(1,2), repel=TRUE)+
  geom_point(aes(Dim.1, Dim.2), data = quali.coord,
             color = "black", size = 1)+ 
  geom_text(aes(Dim.1, Dim.2, label = name), 
            data = quali.coord, vjust = -1, color =  "black", cex=3) + theme_minimal() 

fviz_pca_var(res.pca, 
             invisible=c("ind","var", "quanti.sup"), axes=c(1,3), repel=TRUE)+
  geom_point(aes(Dim.1, Dim.3), data = quali.coord,
             color = "black", size = 1)+ 
  geom_text(aes(Dim.1, Dim.3, label = name), 
            data = quali.coord, vjust = -1, color =  "black", cex=3) + theme_minimal() 

#res.pca.hcpc<-HCPC(res.pca, nb.clust=-1, graph=FALSE) # 6 clusters, but looks strange
#res.pca.hcpc$desc.var
#center_ethic$cluster<-res.pca.hcpc$data.clust[6]

#number in each cluster
#sjt.frq(center_ethic$cluster[6])
#percentage from each country in each cluster


#fviz_cluster(res.pca.hcpc, data = as.numeric(centered_roles), frame.type = "convex", repel=FALSE, geom=c("point"))+ theme_tufte() + labs(title = "C12 clusters (after PCA)") # not working 

#plot.HCPC(res.pca.hcpc, axes=c(1,2), choice="map", ind.names=FALSE, draw.tree = FALSE)
#plot.HCPC(res.pca.hcpc, axes=c(3,4), choice="map", ind.names=FALSE, draw.tree = FALSE)
#plot.HCPC(res.pca.hcpc, axes=c(5,6), choice="map", ind.names=FALSE, draw.tree = FALSE)

#ala SPSS
#pref.res        <- prcomp(ethic[,c(-5)], center=F, scale=T) # returns rotated loadings
#pref.res
#sjt.pca(pref.res, nmbr.fctr = 3, show.msa=TRUE)

#fit <- fa(ethic[,c(-5)], nfactors=4, fm="pa", rotation="varimax")
#fit$loadings
#fit <- fa(center_ethic[,c(-5)], nfactors=4, fm="pa", rotation="varimax")
#fit$loadings

FA test C13

ethic<-BATmC13
#roles<-dplyr::filter(roles,COUNTRY!='Singapore_NON-REP' & COUNTRY !='Ethiopia') # remove Singapore and Ethiopia
# ethic<-ethic %>% filter(complete.cases(.))
x<-t(ethic[,c(-(5:6))])
x<-scale(x, scale=FALSE)
x<-t(x)
center_ethic<-as.data.frame(x)
center_ethic$COUNTRY<-ethic$COUNTRY
center_ethic$ID<-ethic$ID

ethics<-ethic[,-c(5:6)]
#ethics<-center_ethic[,-c(19:20)]

describe(ethics)
KMO(ethics) # not really suitable, KMO MSA=0.67
cortest.bartlett(ethics) # suggest 2 factors and 1 component
fa.parallel(ethics)
scree(ethics)
#varimax, rotated
myfa<-fa(ethics, fm = "pa", nfactors = 2, rotate = "varimax", covar=FALSE) 
print(myfa, cut=.4, sort=TRUE)
plot(myfa, xlim = c(-1, 1), ylim = c(-1, 1))
fa.diagram(myfa, cut = .001, main="FA C13, PCA extraction, varimax")
# oblimin
myfa<-fa(ethics, fm = "pa", nfactors = 2, rotate = "oblimin", covar=FALSE) 
print(myfa, cut=.4, sort=TRUE)
plot(myfa, xlim = c(-1, 1), ylim = c(-1, 1))
fa.diagram(myfa, cut = .001, main="FA C13, PCA extraction, oblimin")



# note: the test does not align exactly with SPSS procedures
#test for scaling
#colab<-dplyr::select(ethics, Pos.image.politicians, Sup.gov.policy)
#alpha(colab)

C14 ETHICS (PRACTICES)

#C14
BATm <- apply(dplyr::select(BAT, Paying.for.info:Accept.money), 2, function(x) {x <- recode(x,"'not approve under any circumstances'=1; 'justified on occasion'=2; 'always justified'=3"); x}) 
### re-add country and add labels
BATm<-data.frame(BATm)
BATm$COUNTRY<-BAT$COUNTRY 
BATm$ID<-BAT$ID 
BATmC14<-BATm
BATmCOMPLETE<-BATm %>% filter(complete.cases(.)) # note filter on only complete C14* cases
BATmORG<-BATm # make a copy of the dataset (orginal)
### convert to matrix of averages by country and add rownames and labels
BATm<-BATm %>%
  group_by(COUNTRY) %>%
  summarise_each(funs(Mean), Paying.for.info:Accept.money)
BATm<-data.frame(BATm)
rownames(BATm)<-BATm[,1]
BATm<-BATm[,2:11]
BATmORGC14<-BATm
### make table of mean and sd by question and country
DT <- data.table(BATmC14)
wide <- setnames(DT[, sapply(.SD, function(x) list(n=sum(complete.cases(x)), mean=round(mean(x, na.rm=TRUE), 2), sd=round(sd(x, na.rm=TRUE), 2))), by=COUNTRY], c("COUNTRY", sapply(names(DT)[-12], paste0, c(".n", ".men", ".SD"))))
table.meansdC14<-as.data.frame(wide[,1:31])

Various descriptive statistics and figures C14

# sjp.setTheme(geom.outline.size = 0, geom.label.size = 3, title.size = 1.5)
sjp.likert(BAT[,6:15], sort.frq = "pos.asc", show.prc.sign = TRUE, digits=0, show.n=FALSE, show.legend=FALSE, catcount=4,
           title="Agreement?     (C14, all countries)") # lickert plot
# Problem: Sweden is missing C14E (Personal documents) and Japan C14J (Publ. unverified)
## solution: The value is set to average for all other countries 
## colMeans(BATm, na.rm=TRUE)
# Japan
BATm[36,9]<-1.236357 #CHECKROWNUMBER 
# Sweden
BATm[60,5]<-1.401936 #CHECKROWNUMBER
# Average mean and spread of roles (based on country averages), boxplot and beanplot
sortedBATm<-BATm[ , order(colMeans(BATm))] # sort by column mean
par(mar = c(4,15,2,2)) # set margins

boxplot(x = as.list(sortedBATm), horizontal = TRUE, las=1, main="Averages aggregated by country (C14)", log="")
par(mar = c(4,15,2,2)) # set margins

beanplot(x = as.list(sortedBATm), horizontal = TRUE, las=1, main="Averages aggregated by country (C14)", log="")

##Bertin plot (centered)
cBATm<-t(BATm)
cBATm<-scale(cBATm, scale=FALSE) #centering of means
cBATm<-t(cBATm)
cBATmC14<-cBATm
x <- scale_by_rank(as.matrix(scale(cBATm, scale=FALSE)))
order <- seriate(x, method="PCA")
bertinplot(x, order, options = list(panel=panel.squares, shading=TRUE, reverse=TRUE, mar = c(1,1,10,6), gp_labels=gpar(fontsize=8)))

order <- seriate(x, method="PCA_angle")
bertinplot(x, order, options = list(panel=panel.squares, shading=TRUE, reverse=TRUE, mar = c(1,1,10,6), gp_labels=gpar(fontsize=8)))

heatmap(cBATm, margins=c(5,5)) # heatmap (centered)

# network plot of correlations
c<-correlate(cBATmC14)
network_plot(c, min_cor=.2, colors=c("red", "green"), legend=TRUE) 

c<-cor(cBATmC14)
corrplot(c, type = "upper", order = "FPC", tl.col = "black", tl.srt = 45, diag=FALSE, tl.cex=0.8,mar=c(1,0,0,1))

#combined plots C13 and C14
c1314<-cbind(cBATmC13,cBATmC14)
c<-cor(c1314)
corrplot(c, type = "upper", order = "original", tl.col = "black", tl.srt = 45, diag=FALSE, tl.cex=0.8,mar=c(1,0,0,1))

c<-correlate(c1314)
network_plot(c, min_cor=.2, colors=c("red", "green"), legend=TRUE)

rplot(shave(c), print_cor=TRUE) + theme(axis.text.x = element_text(angle = 45, hjust = 1))

#canonical correlations
library(CCA)
cc<-cancor(cBATmC13, cBATmC14)
cc

$cor
[1] 0.5991472 0.4495393 0.2199413

$xcoef
                                    [,1]      [,2]         [,3]
Always.follow.codes          -0.15328053 0.4287914  0.002250226
Depends.on.situation          0.01956273 0.5060441  0.550136053
Matter.of.personal.judgement  0.32048287 0.5370145 -0.082250191

$ycoef
                            [,1]        [,2]      [,3]        [,4]          [,5]        [,6]
Paying.for.info       0.12672875  0.61042645 0.8893461  0.69937066  0.0005251803 -0.26478847
Using.conf.documents -0.31936877  0.46739183 0.2615651  0.28401950 -0.2917187026 -0.18598097
False.identity        0.15849511  0.13991719 0.2136405 -0.03739294 -0.3082950217 -0.50111314
Pressure.informants   0.01665925  0.04090378 0.4438876  1.00829851  0.2682033274 -0.25499006
Use.personal.doc      0.72313803 -0.29925869 0.8093177  0.31336474 -1.1992272750 -0.21227944
Undercover.empl       0.17503879 -0.25538778 0.8849925  0.29110893  0.0015249816 -0.80929924
Hidden.miccam        -0.10292762 -0.21966541 0.6530482  0.20240918 -0.0580196397  0.20731753
Recreations           0.07898757  0.10669771 0.7361542  0.19902505 -0.0887596368 -0.02010905
Publ.unverified       0.35536043  0.43101357 0.8706335  0.21687099 -0.0982503683 -0.26274579
                             [,7]         [,8]       [,9]
Paying.for.info       0.319220026 -0.253139609  0.2805626
Using.conf.documents -0.011415615 -0.032758123  0.4183112
False.identity        0.453229994  0.589572658  0.3344843
Pressure.informants   0.119618195  0.039008822  0.9390976
Use.personal.doc      0.017284933  0.368242662 -0.1667630
Undercover.empl      -0.429115367  0.003174688  0.4016016
Hidden.miccam         0.638434324 -0.385980287  0.8567863
Recreations           0.006425329  0.405021116  0.2622131
Publ.unverified       0.269149296 -0.334936280  1.4676279

$xcenter
         Always.follow.codes         Depends.on.situation Matter.of.personal.judgement 
                   1.2367947                   -0.1195097                   -0.5899000 
            Can.be.set.aside 
                  -0.5273851 

$ycenter
     Paying.for.info Using.conf.documents       False.identity  Pressure.informants 
       -0.0233298460         0.2808274257         0.0419988764         0.0003955141 
    Use.personal.doc      Undercover.empl        Hidden.miccam          Recreations 
       -0.1026957262         0.1339708539         0.2616879713         0.0402698942 
     Publ.unverified         Accept.money 
       -0.2690506891        -0.3640742742 

# some nice HTML tables
#sjt.frq(BATm)
#sjt.corr(BATm, triangle="lower", digits=2)

PCA and biplot, mean centered C14

x<-t(cBATm)
x<-scale(x, scale=FALSE)
x<-t(x)
C14x<-x
paran(x)

Using eigendecomposition of correlation matrix.
Computing: 10%  20%  30%  40%  50%  60%  70%  80%  90%  100%


Results of Horn's Parallel Analysis for component retention
300 iterations, using the mean estimate

-------------------------------------------------- 
Component   Adjusted    Unadjusted    Estimated 
            Eigenvalue  Eigenvalue    Bias 
-------------------------------------------------- 
1           2.327468    2.980235      0.652766
2           1.439568    1.877161      0.437593
3           1.173156    1.454701      0.281545
-------------------------------------------------- 

Adjusted eigenvalues > 1 indicate dimensions to retain.
(3 components retained)

res.pca<-PCA(x , scale.unit=FALSE, ncp=3, graph = TRUE)

fviz_pca_biplot(res.pca, geom = "text",  title="Ethics(C14)", axes =c(1,2), labelsize=2.5, repel=TRUE) + theme_minimal()

fviz_pca_biplot(res.pca, geom = "text",  title="Ethics(C14)", axes =c(3,2), labelsize=2.5, repel=TRUE) + theme_minimal()

fviz_screeplot(res.pca, ncp=10)

fviz_contrib(res.pca, choice = "var", axes = 1)

fviz_contrib(res.pca, choice = "var", axes = 2)

fviz_contrib(res.pca, choice = "var", axes = 3)

indroles <- get_pca_ind(res.pca)
C14ind<-indroles$coord
# HCPC Clustering
res.pca.hcpc<-HCPC(res.pca, nb.clust=-1) # 3 clusters

# fviz_cluster(res.pca.hcpc, data = x, frame.type = "convex", repel=TRUE)+ theme_tufte() + labs(title = "C14 clusters (after PCA)") + ylab("2: Paying/being paid --->") + xlab("1: Undercover methods --->") + theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold"))
res.pca.hcpc$desc.var

$quanti.var
                          Eta2      P-value
Using.conf.documents 0.5818686 7.620935e-13
Recreations          0.5787516 9.665527e-13
Accept.money         0.3649914 4.885927e-07
Paying.for.info      0.2375876 1.698551e-04
Publ.unverified      0.2254328 2.817673e-04
False.identity       0.2011706 7.560136e-04
Undercover.empl      0.1657251 3.033151e-03
Hidden.miccam        0.1547437 4.609245e-03

$quanti
$quanti$`1`
                        v.test Mean in category Overall mean sd in category Overall sd
Using.conf.documents  4.346647        0.4126477   0.28082743     0.12169525  0.2152321
Recreations           4.278100        0.1822538   0.04026989     0.15912586  0.2355414
False.identity       -2.668399       -0.0177491   0.04199888     0.13541243  0.1589101
Accept.money         -3.367483       -0.4656382  -0.36407427     0.09507777  0.2140490
Paying.for.info      -3.899232       -0.1105609  -0.02332985     0.11184898  0.1587711
                          p.value
Using.conf.documents 1.382347e-05
Recreations          1.884957e-05
False.identity       7.621358e-03
Accept.money         7.585778e-04
Paying.for.info      9.649843e-05

$quanti$`2`
                   v.test Mean in category Overall mean sd in category Overall sd      p.value
False.identity   3.539530       0.15204355   0.04199888     0.16259073  0.1589101 4.008395e-04
Undercover.empl  2.772809       0.20917751   0.13397085     0.11589843  0.1386323 5.557479e-03
Paying.for.info  2.719943       0.06115974  -0.02332985     0.15237615  0.1587711 6.529313e-03
Hidden.miccam    2.488772       0.34009728   0.26168797     0.17156198  0.1610314 1.281850e-02
Publ.unverified -2.758933      -0.35144996  -0.26905069     0.08896683  0.1526548 5.799044e-03
Recreations     -6.077643      -0.23980511   0.04026989     0.10444094  0.2355414 1.219622e-09

$quanti$`3`
                        v.test Mean in category Overall mean sd in category Overall sd
Accept.money          4.833805      -0.16164463   -0.3640743      0.2822467  0.2140490
Publ.unverified       3.567609      -0.16249921   -0.2690507      0.1803796  0.1526548
Undercover.empl      -2.752943       0.05930301    0.1339709      0.1292516  0.1386323
Hidden.miccam        -2.826162       0.17264913    0.2616880      0.1673959  0.1610314
Using.conf.documents -6.079042       0.02484266    0.2808274      0.1234124  0.2152321
                          p.value
Accept.money         1.339477e-06
Publ.unverified      3.602533e-04
Undercover.empl      5.906214e-03
Hidden.miccam        4.710940e-03
Using.conf.documents 1.209024e-09


attr(,"class")
[1] "catdes" "list " 

# K-means clustering
fviz_nbclust(x, kmeans, method = "gap_stat") # 1 cluster

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 100)  [one "." per sample]:
.................................................. 50 
.................................................. 100 

# Non-metric MDS
library(ggrepel)
d<-dist(C14x)
fit <- isoMDS(d, k=2)

initial  value 26.682526 
iter   5 value 19.082329
iter  10 value 18.047491
final  value 17.916834 
converged

fit

$points
                             [,1]         [,2]
Albania              -0.153980733  0.694066765
Argentina            -0.134522694 -0.085624313
Australia             0.049993565 -0.755301607
Austria              -0.182261063 -0.370090605
Bangladesh            0.048016389  0.397134465
Belgium              -0.361470897 -0.002427200
Bhutan                0.360551581  0.217577221
Botswana             -0.059960785  0.015350845
Brazil               -0.503113321  0.195750654
Bulgaria             -0.350866576  0.226109952
Canada               -0.564892488 -0.166988101
Chile                -0.060972886 -0.024982137
China                 1.089675258  0.288444591
Colombia              0.263718742 -0.202023499
Croatia              -0.085747428 -0.249366798
Cyprus               -0.028583528  0.508923836
Czech Republic       -0.263959025  0.222131613
Denmark              -0.420015903 -0.264637122
Ecuador               0.332570882 -0.349633019
Egypt                 0.113508893  0.422761806
El Salvador           0.123759312 -0.134711292
Estonia               0.157311377 -0.431727567
Ethiopia              0.431248587  0.048515502
Finland              -0.131873310 -0.462744446
France               -0.855428742  0.210186018
Germany              -0.092746701 -0.574885553
Greece               -0.134408613  0.273654345
Hong Kong            -0.006113759 -0.032191944
Hungary               0.230075269  0.084494754
Iceland              -0.444245934  0.048720057
India                 0.142600734 -0.048561276
Indonesia             0.111342241  0.550566669
Ireland              -0.345439396 -0.334793136
Israel               -0.471652972 -0.548260653
Italy                -0.314736417 -0.164902611
Japan                -0.158682132 -0.750891764
Kenya                 0.140373557  0.031652769
Kosovo               -0.040544824  0.239446597
Latvia               -0.309003933  0.193001639
Malawi                0.123829087  0.135735462
Malaysia              0.381634778 -0.025202245
Mexico               -0.080384157 -0.057876794
Moldova              -0.306630231  0.487881535
Netherlands          -0.344174410  0.133183509
New Zealand          -0.107058406 -0.297954950
Norway               -0.327644323 -0.019096215
Oman                  0.724379556 -0.409824945
Philippines           0.018277308 -0.248114003
Portugal              0.286877463 -0.106922789
Qatar                 0.534757180  0.020668407
Romania              -0.217108856  0.538429265
Russia                0.039782183  0.351932567
Serbia                0.015965909  0.030515355
Sierra Leone          0.262577928  0.238593709
Singapore             0.137619159 -0.182261678
South Africa         -0.193798046 -0.105680876
South Korea          -0.226227495  0.375762769
Spain                -0.279244457 -0.013539542
Sudan                 0.191198778  0.468538043
Sweden               -0.471264515 -0.128625322
Switzerland          -0.167019710 -0.324341296
Tanzania              2.999240337 -0.002650837
Thailand              0.448497501 -0.098042357
Turkey               -0.462197812  0.663916482
United Arab Emirates  0.386743856  0.099538650
UK                   -0.106140958 -0.216272049
USA                  -0.382009977 -0.222035310

$stress
[1] 17.91683

x <- fit$points[,1]
y <- fit$points[,2]
plot(x, y, xlab="Coordinate 1", ylab="Coordinate 2",
main="NonMetric MDS",   type="n")
text(x, y, labels = row.names(BATm), cex=.7)

space<-as.data.frame(fit$points)
space$label<-C14x[,-1]
qplot(V1, V2, data=space) + geom_text_repel(aes(label=row.names(fit$points))) + ggtitle("Non-parametric MDS C14") + theme_minimal()

#res.mds.hcpc<-HCPC(space, nb.clust=-1) 
#fviz_cluster(res.mds.hcpc, data = x, frame.type = "convex", repel=TRUE, labelsize=3)+ theme_tufte() + labs(title = "C14 clusters (after MDS)") 
# res.mds.hcpc$desc.var

PCA C14 ETHICS (non-aggregated)

ethic<-BATmC14
#roles<-dplyr::filter(roles,COUNTRY!='Singapore_NON-REP' & COUNTRY !='Ethiopia') # remove Singapore and Ethiopia
ethic<-ethic %>% filter(complete.cases(.))

x<-t(ethic[,-c(11:12)])
x<-scale(x, scale=FALSE)
x<-t(x)
center_ethic<-as.data.frame(x)
center_ethic$COUNTRY<-ethic$COUNTRY
center_ethic$ID<-ethic$ID

KMO(as.matrix(center_ethic[,-c(11:12)]))
paran(as.matrix(center_ethic[,-c(11:12)])) #5 components

# non-rotated solution, with centered data
res.pca<-PCA(center_ethic, quali.sup=c(11:12), scale.unit=TRUE, ncp=5, graph = FALSE) 
summary.PCA(res.pca)
res.pca$var
sweep(res.pca$var$coord,2,sqrt(res.pca$eig[1:ncol(res.pca$var$coord),1]),FUN="/")
center_ethic$dim1<-res.pca$ind[1]
center_ethic$dim2<-res.pca$ind[2]
center_ethic$dim3<-res.pca$ind[3]
center_ethic$dim4<-res.pca$ind[4]
center_ethic$dim5<-res.pca$ind[5]

fviz_screeplot(res.pca, ncp=10)
fviz_contrib(res.pca, choice = "var", axes = 1)
fviz_contrib(res.pca, choice = "var", axes = 2)
fviz_contrib(res.pca, choice = "var", axes = 3)
fviz_contrib(res.pca, choice = "var", axes = 4)
fviz_contrib(res.pca, choice = "var", axes = 5)


res.pca$var

# rotated loadings
loadings<-sweep(res.pca$var$coord,2,sqrt(res.pca$eig[1:5,1]),FUN="/")
loadings
#fviz_pca_biplot(res.pca, geom = "text",  title="Roles(C12)", axes =c(1,2), repel=FALSE) + theme_minimal()
fviz_pca_biplot(res.pca, geom = "text",  invisible="ind", title="Ethics(C14)", axes =c(1,2), repel=TRUE) + theme_minimal()
fviz_pca_biplot(res.pca, geom = "text",  invisible="ind", title="Ethics(C14)", axes =c(3,4), repel=TRUE) + theme_minimal()
fviz_pca_biplot(res.pca, geom = "text",  invisible="ind", title="Ethics(C14)", axes =c(1,5), repel=TRUE) + theme_minimal()


quali.coord <- as.data.frame(res.pca$quali.sup$coord)
quali.coord$name <- rownames(quali.coord)
quali.coord<-transform(quali.coord, Dim.1=Dim.1, Dim.2=Dim.2, Dim.3=Dim.3, Dim.4=Dim.4, Dim.5=Dim.5)
fviz_pca_var(res.pca, 
             invisible=c("ind","var", "quanti.sup"), axes=c(1,2), repel=TRUE)+
  geom_point(aes(Dim.1, Dim.2), data = quali.coord,
             color = "black", size = 1)+ 
  geom_text(aes(Dim.1, Dim.2, label = name), 
            data = quali.coord, vjust = -1, color =  "black", cex=3) + theme_minimal() 

fviz_pca_var(res.pca, 
             invisible=c("ind","var", "quanti.sup"), axes=c(3,4), repel=TRUE)+
  geom_point(aes(Dim.3, Dim.4), data = quali.coord,
             color = "black", size = 1)+ 
  geom_text(aes(Dim.3, Dim.4, label = name), 
            data = quali.coord, vjust = -1, color =  "black", cex=3) + theme_minimal() 

res.pca.hcpc<-HCPC(res.pca, nb.clust=-1, graph=FALSE) # 6 clusters, but looks strange
res.pca.hcpc$desc.var
center_ethic$cluster<-res.pca.hcpc$data.clust


#number in each cluster
sjt.frq(center_ethic$cluster[6])
#percentage from each country in each cluster


#fviz_cluster(res.pca.hcpc, data = as.numeric(centered_roles), frame.type = "convex", repel=FALSE, geom=c("point"))+ theme_tufte() + labs(title = "C12 clusters (after PCA)") # not working 

plot.HCPC(res.pca.hcpc, axes=c(1,2), choice="map", ind.names=FALSE, draw.tree = FALSE)
plot.HCPC(res.pca.hcpc, axes=c(3,4), choice="map", ind.names=FALSE, draw.tree = FALSE)
plot.HCPC(res.pca.hcpc, axes=c(5,6), choice="map", ind.names=FALSE, draw.tree = FALSE)

FA test C14

prethic<-BATmC14

#prethic<-dplyr::filter(roles,COUNTRY!='Singapore_NON-REP' & COUNTRY !='Ethiopia') # remove Singapore and Ethiopia
x<-t(prethic[,c(-(11:12))])
x<-scale(x, scale=FALSE)
x<-t(x)
center_prethic<-as.data.frame(x)
center_prethic$COUNTRY<-ethic$COUNTRY
center_prethic$ID<-prethic$ID

prethics<-prethic[,-c(11:12)] #removing country and ID
prethics<-prethics[,-8] # removing recreations
#prethics<-center_prethic[,-c(11:12)]
#prethics<-prethics %>% filter(complete.cases(.))


describe(prethics)
KMO(prethics) #  suitable, KMO MSA=0.83, .82 without recreations
cortest.bartlett(prethics)
fa.parallel(prethics)  # suggest 4 factors and 2 components (also without recr.)
#varimax, rotated
myfa<-fa(prethics, fm = "pa", nfactors = 4, rotate = "varimax", covar=FALSE, SMC=FALSE) 
print(myfa, cut=.4, sort=TRUE)
plot(myfa, xlim = c(-1, 1), ylim = c(-1, 1))
fa.diagram(myfa, cut = .001, main="FA C14, PCA extraction, varimax")
# oblimin
myfa<-fa(prethics, fm = "pa", nfactors = 4, rotate = "oblimin", covar=FALSE) 
print(myfa, cut=.4, sort=TRUE)
plot(myfa, xlim = c(-1, 1), ylim = c(-1, 1))
fa.diagram(myfa, cut = .001, main="FA C14, PCA extraction, oblimin")

#test for scaling # none are above 0.70, which is problematic

# but Nunnally (1978): When there are a small number of items in the scale (fewer than 10), Cronbach alpha values can be quite small. In this situation it may be better to calculate and report the mean inter-item correlation for the items. Optimal mean inter-item correlation values range from .2 to .4 (as recommended by Briggs & Cheek 1986). 

covert<-dplyr::select(prethics, Undercover.empl, Hidden.miccam, False.identity,Paying.for.info )
alpha(covert) #0.69, average inter-item correlation = average_r=.33 OK

personal<-dplyr::select(prethics, Using.conf.documents, Pressure.informants, Use.personal.doc)
alpha(personal) #0.59, average_r=.33

personal2<-dplyr::select(prethics, Using.conf.documents, Pressure.informants)
alpha(personal2) #0.4

moneyunver<-dplyr::select(prethics, Publ.unverified, Accept.money)
alpha(moneyunver) #0.64, average_r=.48

#measurement invariance # not completed for now
library(lavaan)
prethics<-prethic[,-12] #removing  ID
prethics<-prethics[,-8] # removing recreations
prethics<-prethics %>% filter(complete.cases(.))

HS.model <- '  undercover =~ Undercover.empl + Hidden.miccam + False.identity+Paying.for.info
              moneyunv =~ Accept.money+Publ.unverified
              pressure   =~ Pressure.informants
              docs =~ Using.conf.documents+Use.personal.doc'

summary(fit)
library(semTools)
#measurementInvariance(HS.model, data = prethics, group = "COUNTRY") # not working