WJS Ethics pr dec18, JFH

Install required packages and set working directory

require(foreign)

Warning messages:
1: In scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  EOF within quoted string
2: In scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  EOF within quoted string

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)
require(xlsx)
require(corrplot)
require(ape)
require(grid)
require(matrixStats)
require(ggthemes)
require(data.table)
require(psych)
require(explor)
library(rworldmap)
library(CCA)
library(ggrepel)

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)

`summarise_each()` is deprecated.
Use `summarise_all()`, `summarise_at()` or `summarise_if()` instead.
To map `funs` over a selection of variables, use `summarise_at()`

BATm<-data.frame(BATm)
rownames(BATm)<-BATm[,1]
BATmORGC13map<-BATm
BATm<-BATm[,2:5]
BATmORGC13<-BATm
#BATm<-BATm[c(-61),]# remove Tanzania 
BATm<-BATm[c(-10),]# remove Bulgaria, C13 has been wrongly translated
BATmORGC13<-BATmORGC13[c(-10),]# remove Bulgaria, C13 has been wrongly translated
### 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

Scatterplots C13

fcountry <- merge(BATmORGC13map,country,by="COUNTRY")
rownames(fcountry)<-fcountry[,1]
namn<-row.names(fcountry) 
fcountry<-fcountry[c(-10),] #removing Bulgaria because of C13 wrong translation
fcountry<-fcountry[c(-22),] #removing Ethiopia = outlier
fcountry<-fcountry[c(-53),] #removing Singapore = outlier
qplot(FHPFREES,EIUDEMOS, data = fcountry) + geom_text_repel(aes(label=row.names(fcountry))) +  theme_minimal()

qplot(FHPFREES,fcountry[,2], data = fcountry) + geom_text_repel(aes(label=row.names(fcountry))) +  theme_minimal() + labs(x = "Freedom House Press Freedom index", y="Always follow codes")

qplot(FHPFREES,fcountry[,3], data = fcountry) + geom_text_repel(aes(label=row.names(fcountry))) +  theme_minimal() + labs(x = "Freedom House Press Freedom index", y="Depends on situation")

qplot(FHPFREES,fcountry[,4], data = fcountry) + geom_text_repel(aes(label=row.names(fcountry))) +  theme_minimal() + labs(x = "Freedom House Press Freedom index", y="Matter of personal judgement")

qplot(FHPFREES,fcountry[,5], data = fcountry) + geom_text_repel(aes(label=row.names(fcountry))) +  theme_minimal() + labs(x = "Freedom House Press Freedom index", y="Can be set aside")

qplot(EIUDEMOS,fcountry[,2], data = fcountry) + geom_text_repel(aes(label=row.names(fcountry))) +  theme_minimal() + labs(x = "EIU Democracy Score", y="Always follow codes")

qplot(EIUDEMOS,fcountry[,3], data = fcountry) + geom_text_repel(aes(label=row.names(fcountry))) +  theme_minimal() + labs(x = "EIU Democracy Score", y="Depends on situation")

qplot(EIUDEMOS,fcountry[,4], data = fcountry) + geom_text_repel(aes(label=row.names(fcountry))) +  theme_minimal() + labs(x = "EIU Democracy Score", y="Matter of personal judgement")

qplot(EIUDEMOS,fcountry[,5], data = fcountry) + geom_text_repel(aes(label=row.names(fcountry))) +  theme_minimal() + labs(x = "EIU Democracy Score", y="Can be set aside")

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)))

#uncentered C13 Bertin
x <- scale_by_rank(as.matrix(BATm))
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)))

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))

sjp.corr(c)

Computing correlation using pearson-method with listwise-deletion...

# distance matrix (correlations)
#res.dist <- get_dist(C13x, stand = TRUE, method = "pearson")
#fviz_dist(res.dist, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))

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
# test for PCA factors = 2
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.611300    1.903982      0.292681
2           1.151001    1.218220      0.067218
-------------------------------------------------- 

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

res.pca<-PCA(x , scale.unit=FALSE, ncp=3, graph = TRUE) # suggests 2 axes only

#res.pca$var$coord
fviz_pca_var(res.pca, col.var = "black")+ labs(title = "C13 PCA, centered means")

#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_pca_biplot(res.pca, repel = TRUE)+ labs(title = "C13 PCA, centered means")

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) # not stable (PARAN)

indroles <- get_pca_ind(res.pca)
C13ind<-indroles$coord
# Quality of representation of the Countries
fviz_cos2(res.pca, choice="ind", axes = 1, top = 68 )

fviz_cos2(res.pca, choice="ind", axes = 2, top = 68 )

# K-menas clustering = just one cluser
set.seed(123)
fviz_nbclust(C13x, 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 

# km.res <- kmeans((C13x), 4, nstart = 25)
#fviz_cluster(km.res, C13x, ellipse.type = "norm")
# HCPC Clustering
res.pca.hcpc<-HCPC(res.pca, nb.clust=-1) # 5 clusters

res.pca.hcpc$data.clust

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

#fviz_cluster(res.pca.hcpc, data = x, axes=c(1,3), ellipse.type = "convex", repel=TRUE)+ theme_tufte() + labs(title = "C13 clusters (after PCA)") + ylab(" <--- 1: Absolutism") + xlab("2: <--- Situationism") + theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold"))
#fviz_cluster(res.pca.hcpc, data = x, axes=c(2,3), ellipse.type = "convex", repel=TRUE)+ theme_tufte() + labs(title = "C13 clusters (after PCA)") + ylab(" <--- 2: Subjectivism") + xlab("3: <--- Situationism") + theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold"))
fviz_dend(res.pca.hcpc, rect = TRUE, cex = 0.5, horiz=TRUE)

#res.pca.hcpc$desc.var
# Non-metric MDS
library(ggrepel)
d<-dist(C13x[,c(-2)]) # also removing Argentina #CHECKROWNUMBER
fit <- isoMDS(d, k=2)

initial  value 4.168022 
iter   5 value 3.861572
final  value 3.850813 
converged

fit

$points
                             [,1]         [,2]
Albania               0.254043740  0.013908584
Argentina            -0.190321897 -0.206896483
Australia            -0.391511646  0.192218283
Austria              -0.629906667 -0.181398664
Bangladesh           -0.211676875 -0.427445210
Belgium              -0.014701740  0.244288859
Bhutan                0.396023651 -0.027502323
Botswana             -0.356831124  0.387599788
Brazil               -0.273172084  0.060477270
Canada               -0.185745158  0.360944338
Chile                -0.265965777  0.028929172
China                -0.218128203  0.155833286
Colombia             -0.363921451 -0.517614539
Croatia              -0.508344660 -0.358308731
Cyprus               -0.262750459  0.434399371
Czech Republic        0.496559686  0.534575410
Denmark               0.536104967  0.350898749
Ecuador              -0.212861654 -0.350732579
Egypt                 0.072108277 -0.566981270
El Salvador          -0.227894360 -0.191658034
Estonia              -0.325734053  0.055164819
Ethiopia              1.501558558  0.440744278
Finland              -0.584368019 -0.278971484
France                0.142208639  0.299972366
Germany              -0.641366091  0.007237585
Greece               -0.143118737  0.287963860
Hong Kong             0.306621330 -0.172937977
Hungary               0.189998006  0.230587975
Iceland              -0.097081313 -0.009127974
India                 0.599173572 -0.033770224
Indonesia            -0.504458948  0.191826206
Ireland               0.118064473  0.049483701
Israel                0.061751627  0.247962909
Italy                -0.561181871  0.343091207
Japan                 0.378287542  0.045154544
Kenya                 0.054513631  0.150248968
Kosovo               -0.282838647  0.087509762
Latvia               -0.067263776 -0.157906085
Malawi               -0.091431050  0.296126245
Malaysia              0.422012457 -0.211374453
Mexico               -0.002635776 -0.093035264
Moldova              -0.154476411  0.102030449
Netherlands           0.461258492  0.006643104
New Zealand          -0.182913021 -0.012346573
Norway               -0.186477321 -0.227743120
Oman                  0.822149720 -0.491811532
Philippines           0.035349359 -0.191440793
Portugal             -0.477727144  0.315342435
Qatar                 0.373255946 -0.597499952
Romania              -0.298003863  0.092976317
Russia                0.459319604 -0.054645407
Serbia               -0.618596061 -0.312718765
Sierra Leone          0.326634824  0.444646596
Singapore             1.328001454  0.020847490
South Africa         -0.253609451  0.078238161
South Korea          -0.305120500 -0.145853967
Spain                -0.366342896  0.016166388
Sudan                 0.764508431 -0.637250072
Sweden                0.392414266  0.350079042
Switzerland          -0.343309475  0.156814509
Tanzania              0.631313795 -0.149594921
Thailand              0.728207570 -0.162978202
Turkey               -0.321923059 -0.169862138
United Arab Emirates -0.077517746 -0.330048269
UK                    0.012584327  0.005249021
USA                  -0.662798961  0.183273960

$stress
[1] 3.850813

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, centered means") + 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

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)

                                                                  25%           50%           75%                             
 1.079350e+05  3.209450e+00  1.464182e+00  1.000000e+00  2.000000e+00  4.000000e+00  4.000000e+00  5.000000e+00 -1.079310e+05 

KMO(ethics) # not really suitable, KMO MSA=0.67

Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = ethics)
Overall MSA =  0.67
MSA for each item = 
         Always.follow.codes         Depends.on.situation Matter.of.personal.judgement             Can.be.set.aside 
                        0.71                         0.65                         0.66                         0.70 

cortest.bartlett(ethics) # suggest 2 factors and 1 component

R was not square, finding R from data

$chisq
[1] 14166.94

$p.value
[1] 0

$df
[1] 6

fa.parallel(ethics)

Parallel analysis suggests that the number of factors =  2  and the number of components =  1 

scree(ethics)

#varimax, rotated
myfa<-fa(ethics, fm = "pa", nfactors = 2, rotate = "varimax", covar=FALSE) 
print(myfa, cut=.4, sort=TRUE)

Factor Analysis using method =  pa
Call: fa(r = ethics, nfactors = 2, rotate = "varimax", covar = FALSE, 
    fm = "pa")
Standardized loadings (pattern matrix) based upon correlation matrix
                             item   PA1   PA2   h2   u2 com
Matter.of.personal.judgement    3  0.67       0.45 0.55 1.0
Depends.on.situation            2  0.64       0.49 0.51 1.4
Can.be.set.aside                4  0.46  0.42 0.39 0.61 2.0
Always.follow.codes             1       -0.42 0.19 0.81 1.1

                       PA1  PA2
SS loadings           1.07 0.45
Proportion Var        0.27 0.11
Cumulative Var        0.27 0.38
Proportion Explained  0.71 0.29
Cumulative Proportion 0.71 1.00

Mean item complexity =  1.4
Test of the hypothesis that 2 factors are sufficient.

The degrees of freedom for the null model are  6  and the objective function was  0.51 with Chi Square of  14166.94
The degrees of freedom for the model are -1  and the objective function was  0 

The root mean square of the residuals (RMSR) is  0 
The df corrected root mean square of the residuals is  NA 

The harmonic number of observations is  26803 with the empirical chi square  0.44  with prob <  NA 
The total number of observations was  27567  with Likelihood Chi Square =  0.58  with prob <  NA 

Tucker Lewis Index of factoring reliability =  1.001
Fit based upon off diagonal values = 1
Measures of factor score adequacy             
                                                   PA1   PA2
Correlation of (regression) scores with factors   0.78  0.56
Multiple R square of scores with factors          0.61  0.31
Minimum correlation of possible factor scores     0.22 -0.38

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) 

Loading required namespace: GPArotation

print(myfa, cut=.4, sort=TRUE)

Factor Analysis using method =  pa
Call: fa(r = ethics, nfactors = 2, rotate = "oblimin", covar = FALSE, 
    fm = "pa")
Standardized loadings (pattern matrix) based upon correlation matrix
                             item  PA1   PA2   h2   u2 com
Matter.of.personal.judgement    3 0.72       0.45 0.55 1.0
Depends.on.situation            2 0.64       0.49 0.51 1.1
Can.be.set.aside                4 0.41       0.39 0.61 1.9
Always.follow.codes             1       0.43 0.19 0.81 1.0

                       PA1  PA2
SS loadings           1.15 0.37
Proportion Var        0.29 0.09
Cumulative Var        0.29 0.38
Proportion Explained  0.75 0.25
Cumulative Proportion 0.75 1.00

 With factor correlations of 
      PA1   PA2
PA1  1.00 -0.47
PA2 -0.47  1.00

Mean item complexity =  1.3
Test of the hypothesis that 2 factors are sufficient.

The degrees of freedom for the null model are  6  and the objective function was  0.51 with Chi Square of  14166.94
The degrees of freedom for the model are -1  and the objective function was  0 

The root mean square of the residuals (RMSR) is  0 
The df corrected root mean square of the residuals is  NA 

The harmonic number of observations is  26803 with the empirical chi square  0.44  with prob <  NA 
The total number of observations was  27567  with Likelihood Chi Square =  0.58  with prob <  NA 

Tucker Lewis Index of factoring reliability =  1.001
Fit based upon off diagonal values = 1
Measures of factor score adequacy             
                                                   PA1   PA2
Correlation of (regression) scores with factors   0.83  0.64
Multiple R square of scores with factors          0.69  0.40
Minimum correlation of possible factor scores     0.38 -0.19

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)

`summarise_each()` is deprecated.
Use `summarise_all()`, `summarise_at()` or `summarise_if()` instead.
To map `funs` over a selection of variables, use `summarise_at()`

BATm<-data.frame(BATm)
rownames(BATm)<-BATm[,1]
BATm<-BATm[,2:11]
BATm<-BATm[c(-36),] #removing Japan because of missing C14 question
BATm<-BATm[c(-59),] #removing Sweden because of missing C14 question
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])

PCA and biplot, raw values C14

x<-BATm
paran(x) # 2 components

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           3.590664    4.250436      0.659772
2           1.400674    1.839303      0.438628
-------------------------------------------------- 

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

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

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

fviz_pca_var(res.pca, col.var = "black")+ labs(title = "Ethics(C14) Raw values")

fviz_screeplot(res.pca, ncp=10)

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

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

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: Accept money/publ. unverified --->") + xlab("1: Acceptance (general) --->") + theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold"))

argument frame is deprecated; please use ellipse instead.argument frame.type is deprecated; please use ellipse.type instead.

res.pca.hcpc$desc.var

$quanti.var
                          Eta2      P-value
False.identity       0.5054515 3.314947e-10
Undercover.empl      0.4514724 8.225306e-09
Hidden.miccam        0.4499298 8.973496e-09
Using.conf.documents 0.3914928 2.052281e-07
Publ.unverified      0.3874586 2.518830e-07
Recreations          0.3493974 1.632261e-06
Use.personal.doc     0.3323994 3.630936e-06
Pressure.informants  0.2894590 2.507629e-05
Paying.for.info      0.2814599 3.547986e-05
Accept.money         0.2220268 4.167763e-04

$quanti
$quanti$`1`
                        v.test Mean in category Overall mean sd in category Overall sd      p.value
Pressure.informants  -2.056833         1.433895     1.504220      0.1684411  0.1956503 3.970233e-02
Publ.unverified      -2.397120         1.151667     1.230313      0.1018828  0.1877376 1.652450e-02
Paying.for.info      -3.074533         1.369469     1.478706      0.1649810  0.2033088 2.108326e-03
Recreations          -3.464035         1.355207     1.530467      0.2603903  0.2895131 5.321373e-04
Use.personal.doc     -3.659470         1.301948     1.399896      0.1206191  0.1531598 2.527375e-04
Using.conf.documents -4.787314         1.557939     1.773167      0.2461223  0.2572608 1.690284e-06
False.identity       -5.181118         1.338812     1.544273      0.1427306  0.2269199 2.205596e-07
Undercover.empl      -5.206467         1.443252     1.633863      0.1490712  0.2094938 1.924699e-07
Hidden.miccam        -5.335515         1.542969     1.760803      0.2013912  0.2336225 9.527376e-08

$quanti$`2`
                        v.test Mean in category Overall mean sd in category Overall sd      p.value
Using.conf.documents  4.265682         1.921332     1.773167     0.16501165  0.2572608 1.992924e-05
Hidden.miccam         3.104298         1.858720     1.760803     0.14528254  0.2336225 1.907314e-03
Undercover.empl       2.446191         1.703053     1.633863     0.16151870  0.2094938 1.443746e-02
Accept.money         -2.798442         1.064234     1.141770     0.07363577  0.2052133 5.134973e-03

$quanti$`3`
                      v.test Mean in category Overall mean sd in category Overall sd      p.value
Publ.unverified     4.933383         1.461858     1.230313      0.2441595  0.1877376 8.081738e-07
Pressure.informants 4.266388         1.712900     1.504220      0.1839708  0.1956503 1.986635e-05
Recreations         4.244964         1.837710     1.530467      0.2274297  0.2895131 2.186288e-05
False.identity      4.044050         1.773692     1.544273      0.1255413  0.2269199 5.253563e-05
Use.personal.doc    3.930550         1.550396     1.399896      0.1768360  0.1531598 8.475184e-05
Paying.for.info     3.831924         1.673472     1.478706      0.2263899  0.2033088 1.271453e-04
Accept.money        3.534094         1.323081     1.141770      0.2173916  0.2052133 4.091754e-04
Undercover.empl     3.110394         1.796765     1.633863      0.1501910  0.2094938 1.868377e-03
Hidden.miccam       2.442858         1.903479     1.760803      0.1810163  0.2336225 1.457148e-02


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(BATm) # also removing Argentina #CHECKROWNUMBER
fit <- isoMDS(d, k=2)

initial  value 21.272747 
iter   5 value 15.737199
iter  10 value 15.148189
iter  15 value 14.944642
iter  15 value 14.934252
iter  15 value 14.934252
final  value 14.934252 
converged

fit

$points
                             [,1]         [,2]
Albania               0.032773678 -0.518224621
Argentina            -0.184649319 -0.013511940
Australia            -0.275344982  0.483536939
Austria              -0.002365233  0.168018257
Bangladesh           -0.540787403 -0.257746004
Belgium               0.201331779 -0.143793225
Bhutan                0.800739794  0.192704727
Botswana             -0.153520101 -0.055115182
Brazil               -0.103242458 -0.368165731
Bulgaria             -0.004116274 -0.264569501
Canada                0.368404062 -0.208664097
Chile                -0.005983910 -0.033920615
China                 0.489474848  0.835555187
Colombia             -0.711095082  0.191822488
Croatia               0.047797534  0.108327914
Cyprus               -0.246612834 -0.349204182
Czech Republic        0.280161229 -0.228679416
Denmark               0.521694586 -0.021940383
Ecuador              -0.322881918  0.345268862
Egypt                -0.299711681 -0.281387772
El Salvador          -0.455797055  0.110180560
Estonia               0.352799003  0.317314665
Ethiopia             -0.611387880  0.115513122
Finland               0.404591713  0.240081496
France                0.242004486 -0.594867767
Germany               0.152427374  0.357355884
Greece               -0.397536703 -0.190011799
Hong Kong             0.333219194  0.051526775
Hungary               0.258394886  0.170249628
Iceland              -0.155680513 -0.201158078
India                 0.078564284  0.089396777
Indonesia            -0.362256295 -0.405859676
Ireland               0.442862885 -0.000684754
Israel                1.370728591  0.002644432
Italy                 0.232798874 -0.030263571
Kenya                 0.111231202  0.065399486
Kosovo               -0.211302563 -0.189143077
Latvia                0.275051856 -0.200795797
Malawi               -0.004653981 -0.013352845
Malaysia             -0.250756105  0.212361309
Mexico               -0.179652724  0.010582882
Moldova               0.285879350 -0.401862620
Netherlands           0.222978382 -0.221224082
New Zealand          -0.083339715  0.126427159
Norway                0.219132742 -0.120947999
Oman                  1.116528311  0.782749910
Philippines          -0.340005426  0.126085042
Portugal              1.017167122  0.355618988
Qatar                -1.440800825 -0.032913338
Romania               0.327171867 -0.394275466
Russia                0.546988540 -0.181193406
Serbia               -0.231377697 -0.014675825
Sierra Leone          0.471922802  0.144307796
Singapore            -0.410216149  0.135892365
South Africa         -0.116228145 -0.036516379
South Korea           0.112004158 -0.302225828
Spain                -0.047865862 -0.147373888
Sudan                -0.447171933 -0.287793182
Switzerland          -0.031294811  0.171696122
Tanzania             -1.706010536  0.937779542
Thailand              0.161719293  0.416745065
Turkey               -0.018553723 -0.603575024
United Arab Emirates -0.881743218  0.019871498
UK                    0.155482949  0.087514680
USA                  -0.400084319 -0.056892489

$stress
[1] 14.93425

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 C14, centered means") + 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

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 (uncentered)
x <- scale_by_rank(as.matrix(BATm))
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)))

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

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))

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)

                                                                  25%           50%           75%                             
 2.374840e+05  1.501116e+00  6.148928e-01  1.000000e+00  1.000000e+00  1.000000e+00  2.000000e+00  3.000000e+00 -2.374750e+05 

KMO(prethics) #  suitable, KMO MSA=0.83, .82 without recreations

Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = prethics)
Overall MSA =  0.82
MSA for each item = 
     Paying.for.info Using.conf.documents       False.identity  Pressure.informants     Use.personal.doc      Undercover.empl 
                0.87                 0.82                 0.86                 0.88                 0.85                 0.84 
       Hidden.miccam      Publ.unverified         Accept.money 
                0.83                 0.73                 0.68 

cortest.bartlett(prethics)

R was not square, finding R from data

$chisq
[1] 47616.4

$p.value
[1] 0

$df
[1] 36

fa.parallel(prethics)  # suggest 4 factors and 2 components (also without recr.)

Parallel analysis suggests that the number of factors =  4  and the number of components =  2 

#varimax, rotated
myfa<-fa(prethics, fm = "pa", nfactors = 4, rotate = "varimax", covar=FALSE, SMC=FALSE) 
print(myfa, cut=.4, sort=TRUE)

Factor Analysis using method =  pa
Call: fa(r = prethics, nfactors = 4, rotate = "varimax", SMC = FALSE, 
    covar = FALSE, fm = "pa")
Standardized loadings (pattern matrix) based upon correlation matrix
                     item  PA1  PA2  PA3   PA4   h2   u2 com
Hidden.miccam           7 0.62                 0.43 0.57 1.2
Undercover.empl         6 0.62                 0.45 0.55 1.4
False.identity          3 0.57                 0.40 0.60 1.5
Paying.for.info         1 0.41                 0.22 0.78 1.6
Accept.money            9      0.80            0.67 0.33 1.1
Publ.unverified         8      0.57            0.38 0.62 1.4
Pressure.informants     4           0.88       0.87 0.13 1.2
Using.conf.documents    2                 0.54 0.42 0.58 1.8
Use.personal.doc        5                 0.43 0.40 0.60 3.2

                       PA1  PA2  PA3  PA4
SS loadings           1.54 1.16 0.88 0.67
Proportion Var        0.17 0.13 0.10 0.07
Cumulative Var        0.17 0.30 0.40 0.47
Proportion Explained  0.36 0.27 0.21 0.16
Cumulative Proportion 0.36 0.63 0.84 1.00

Mean item complexity =  1.6
Test of the hypothesis that 4 factors are sufficient.

The degrees of freedom for the null model are  36  and the objective function was  1.73 with Chi Square of  47616.4
The degrees of freedom for the model are 6  and the objective function was  0.01 

The root mean square of the residuals (RMSR) is  0.01 
The df corrected root mean square of the residuals is  0.02 

The harmonic number of observations is  25883 with the empirical chi square  176.34  with prob <  2e-35 
The total number of observations was  27567  with Likelihood Chi Square =  239.84  with prob <  6.1e-49 

Tucker Lewis Index of factoring reliability =  0.971
RMSEA index =  0.038  and the 90 % confidence intervals are  0.034 0.042
BIC =  178.5
Fit based upon off diagonal values = 1
Measures of factor score adequacy             
                                                   PA1  PA2  PA3   PA4
Correlation of (regression) scores with factors   0.79 0.84 0.91  0.63
Multiple R square of scores with factors          0.62 0.70 0.82  0.40
Minimum correlation of possible factor scores     0.24 0.41 0.64 -0.20

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,  SMC=FALSE) 
print(myfa, cut=.4, sort=TRUE)

Factor Analysis using method =  pa
Call: fa(r = prethics, nfactors = 4, rotate = "oblimin", SMC = FALSE, 
    covar = FALSE, fm = "pa")
Standardized loadings (pattern matrix) based upon correlation matrix
                     item   PA1   PA2   PA3   PA4   h2   u2 com
Hidden.miccam           7  0.66                   0.43 0.57 1.0
Undercover.empl         6  0.64                   0.45 0.55 1.0
False.identity          3  0.57                   0.40 0.60 1.1
Paying.for.info         1  0.43                   0.22 0.78 1.3
Accept.money            9        0.82             0.67 0.33 1.0
Publ.unverified         8        0.55             0.38 0.62 1.4
Pressure.informants     4              0.93       0.87 0.13 1.0
Using.conf.documents    2                    0.53 0.42 0.58 1.3
Use.personal.doc        5                    0.42 0.40 0.60 1.8

                       PA1  PA2  PA3  PA4
SS loadings           1.53 1.11 0.95 0.66
Proportion Var        0.17 0.12 0.11 0.07
Cumulative Var        0.17 0.29 0.40 0.47
Proportion Explained  0.36 0.26 0.22 0.15
Cumulative Proportion 0.36 0.62 0.85 1.00

 With factor correlations of 
     PA1  PA2  PA3  PA4
PA1 1.00 0.29 0.42 0.61
PA2 0.29 1.00 0.29 0.19
PA3 0.42 0.29 1.00 0.37
PA4 0.61 0.19 0.37 1.00

Mean item complexity =  1.2
Test of the hypothesis that 4 factors are sufficient.

The degrees of freedom for the null model are  36  and the objective function was  1.73 with Chi Square of  47616.4
The degrees of freedom for the model are 6  and the objective function was  0.01 

The root mean square of the residuals (RMSR) is  0.01 
The df corrected root mean square of the residuals is  0.02 

The harmonic number of observations is  25883 with the empirical chi square  176.34  with prob <  2e-35 
The total number of observations was  27567  with Likelihood Chi Square =  239.84  with prob <  6.1e-49 

Tucker Lewis Index of factoring reliability =  0.971
RMSEA index =  0.038  and the 90 % confidence intervals are  0.034 0.042
BIC =  178.5
Fit based upon off diagonal values = 1
Measures of factor score adequacy             
                                                   PA1  PA2  PA3  PA4
Correlation of (regression) scores with factors   0.86 0.86 0.93 0.77
Multiple R square of scores with factors          0.75 0.74 0.87 0.59
Minimum correlation of possible factor scores     0.49 0.47 0.74 0.17

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

Reliability analysis   
Call: alpha(x = covert)

  raw_alpha std.alpha G6(smc) average_r S/N    ase mean   sd
      0.69      0.69    0.63      0.36 2.2 0.0031  1.6 0.44

 lower alpha upper     95% confidence boundaries
0.68 0.69 0.69 

 Reliability if an item is dropped:
                raw_alpha std.alpha G6(smc) average_r S/N alpha se
Undercover.empl      0.59      0.59    0.50      0.33 1.5   0.0042
Hidden.miccam        0.61      0.61    0.51      0.34 1.5   0.0041
False.identity       0.60      0.60    0.51      0.33 1.5   0.0042
Paying.for.info      0.69      0.69    0.59      0.42 2.2   0.0033

 Item statistics 
                    n raw.r std.r r.cor r.drop mean   sd
Undercover.empl 26361  0.76  0.75  0.63   0.52  1.7 0.62
Hidden.miccam   26655  0.73  0.74  0.61   0.50  1.8 0.59
False.identity  26634  0.74  0.74  0.61   0.51  1.5 0.60
Paying.for.info 26210  0.65  0.65  0.44   0.37  1.5 0.60

Non missing response frequency for each item
                   1    2    3 miss
Undercover.empl 0.43 0.49 0.08 0.04
Hidden.miccam   0.31 0.60 0.09 0.03
False.identity  0.51 0.43 0.06 0.03
Paying.for.info 0.58 0.36 0.05 0.05

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

Reliability analysis   
Call: alpha(x = personal)

  raw_alpha std.alpha G6(smc) average_r S/N    ase mean   sd
      0.58      0.59    0.49      0.32 1.4 0.0044  1.6 0.45

 lower alpha upper     95% confidence boundaries
0.57 0.58 0.59 

 Reliability if an item is dropped:
                     raw_alpha std.alpha G6(smc) average_r  S/N alpha se
Using.conf.documents      0.52      0.53    0.36      0.36 1.11   0.0057
Pressure.informants       0.53      0.53    0.36      0.36 1.13   0.0056
Use.personal.doc          0.39      0.39    0.24      0.24 0.64   0.0073

 Item statistics 
                         n raw.r std.r r.cor r.drop mean   sd
Using.conf.documents 26532  0.74  0.72  0.48   0.36  1.8 0.64
Pressure.informants  26441  0.73  0.72  0.48   0.36  1.5 0.62
Use.personal.doc     26013  0.75  0.77  0.60   0.46  1.4 0.56

Non missing response frequency for each item
                        1    2    3 miss
Using.conf.documents 0.32 0.55 0.13 0.04
Pressure.informants  0.57 0.37 0.06 0.04
Use.personal.doc     0.64 0.32 0.04 0.06

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

Reliability analysis   
Call: alpha(x = personal2)

  raw_alpha std.alpha G6(smc) average_r  S/N    ase mean  sd
      0.39      0.39    0.24      0.24 0.64 0.0073  1.7 0.5

 lower alpha upper     95% confidence boundaries
0.38 0.39 0.4 

 Reliability if an item is dropped:
                     raw_alpha std.alpha G6(smc) average_r S/N alpha se
Using.conf.documents      0.24      0.24   0.059      0.24  NA       NA
Pressure.informants       0.24      0.24   0.059      0.24  NA       NA

 Item statistics 
                         n raw.r std.r r.cor r.drop mean   sd
Using.conf.documents 26532  0.80  0.79  0.39   0.24  1.8 0.64
Pressure.informants  26441  0.78  0.79  0.39   0.24  1.5 0.62

Non missing response frequency for each item
                        1    2    3 miss
Using.conf.documents 0.32 0.55 0.13 0.04
Pressure.informants  0.57 0.37 0.06 0.04

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

Reliability analysis   
Call: alpha(x = moneyunver)

  raw_alpha std.alpha G6(smc) average_r S/N    ase mean   sd
      0.63      0.64    0.47      0.47 1.8 0.0043  1.2 0.39

 lower alpha upper     95% confidence boundaries
0.62 0.63 0.64 

 Reliability if an item is dropped:
                raw_alpha std.alpha G6(smc) average_r S/N alpha se
Publ.unverified      0.47      0.47    0.22      0.47  NA       NA
Accept.money         0.47      0.47    0.22      0.47  NA       NA

 Item statistics 
                    n raw.r std.r r.cor r.drop mean   sd
Publ.unverified 25920  0.89  0.86  0.59   0.47  1.2 0.50
Accept.money    26718  0.83  0.86  0.59   0.47  1.1 0.41

Non missing response frequency for each item
                   1    2    3 miss
Publ.unverified 0.79 0.18 0.03 0.06
Accept.money    0.92 0.05 0.03 0.03