#install packages
install.packages(c("psych","car","flextable","officer","data.table","dplyr","numform","effectsize","ppcor"))


#load packages
library(psych)
library(car)
library(flextable)
library(officer)
library(data.table)
library(dplyr)
library(numform)
library(effectsize)
library(ppcor)


#save directory to the folder to which files created in R should be saved (e.g., "C:/NarcSub/")
files_wd <- "C:/NarcSub/"


#read data
dat1 <- as.data.frame(fread("https://madata.bib.uni-mannheim.de/427/3/NarcSub_Study1_Data.csv", header = T, sep = ","))


#apply data-exclusion criteria
dat1 <- subset(dat1, subset = participated_seriously != 2 & not_use_data != 1)


#sample demographics
##sex
dat1_1 <- table(dat1$sex)
length(which(is.na(dat1$sex))) #no missing values
##print sample demographics
paste0(dat1_1[2], " women, ", dat1_1[1], " men; age: ", min(dat1$age), "-", max(dat1$age)," years, M = " , (round(mean(dat1$age),2)), ", SD = ", (round(sd(dat1$age),2)))
rm(dat1_1)


#recode variables
dat1$im1r <- car::recode(dat1$im1, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")
dat1$im3r <- car::recode(dat1$im3, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")
dat1$im4r <- car::recode(dat1$im4, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")
dat1$im5r <- car::recode(dat1$im5, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")
dat1$description_aut_r <- car::recode(dat1$description_aut, "1=6; 2=5; 3=4; 4=3; 5=2; 6=1")
dat1$lbdq2r <- car::recode(dat1$lbdq2, "1=6; 2=5; 3=4; 4=3; 5=2; 6=1")
dat1$lbdq4r <- car::recode(dat1$lbdq4, "1=6; 2=5; 3=4; 4=3; 5=2; 6=1")
dat1$lbdq6r <- car::recode(dat1$lbdq6, "1=6; 2=5; 3=4; 4=3; 5=2; 6=1")
dat1$lbdq8r <- car::recode(dat1$lbdq8, "1=6; 2=5; 3=4; 4=3; 5=2; 6=1")
dat1$lbdq10r <- car::recode(dat1$lbdq10, "1=6; 2=5; 3=4; 4=3; 5=2; 6=1")
dat1$lbdq12r <- car::recode(dat1$lbdq12, "1=6; 2=5; 3=4; 4=3; 5=2; 6=1")
dat1$lbdq14r <- car::recode(dat1$lbdq14, "1=6; 2=5; 3=4; 4=3; 5=2; 6=1")
dat1$lbdq16r <- car::recode(dat1$lbdq16, "1=6; 2=5; 3=4; 4=3; 5=2; 6=1")
dat1$lbdq18r <- car::recode(dat1$lbdq18, "1=6; 2=5; 3=4; 4=3; 5=2; 6=1")
dat1$lbdq20r <- car::recode(dat1$lbdq20, "1=6; 2=5; 3=4; 4=3; 5=2; 6=1")


#compute variables
dat1$narq_adm <- rowMeans( dat1[c("narq1", "narq2", "narq3", "narq5", "narq7", "narq8", "narq15", "narq16", "narq18")], na.rm=TRUE)
dat1$Znarq_adm <- scale(dat1$narq_adm)
dat1$narq_riv <- rowMeans( dat1[c("narq4", "narq6", "narq9", "narq10", "narq11", "narq12", "narq13", "narq14", "narq17")], na.rm=TRUE)
dat1$Znarq_riv <- scale(dat1$narq_riv)
dat1$IM <- rowMeans( dat1[c("im1r","im2","im3r","im4r","im5r","im6","im7","im8")], na.rm=TRUE)
dat1$ZIM <- scale(dat1$IM)
dat1$descriptions <- rowMeans(dat1[c("description_dem","description_aut_r")], na.rm=TRUE)
dat1$Zdescriptions <- scale(dat1$descriptions)
dat1$lbdq <- rowMeans(dat1[c("lbdq1","lbdq3","lbdq5","lbdq7","lbdq9","lbdq11","lbdq13","lbdq15","lbdq17","lbdq19", 
                             "lbdq2r","lbdq4r","lbdq6r","lbdq8r","lbdq10r","lbdq12r","lbdq14r","lbdq16r","lbdq18r","lbdq20r")], na.rm=TRUE)
dat1$Zlbdq <- scale(dat1$lbdq)
dat1$pref_index <- rowMeans(dat1[c("Zdescriptions","Zlbdq")], na.rm=TRUE)
dat1$Zpref_index <- scale(dat1$pref_index)


#internal consistencies of variables
dat.narq_adm <- subset(dat1, select = c(narq1,narq2,narq3,narq5,narq7,narq8,narq15,narq16,narq18))
dat.narq_riv <- subset(dat1, select = c(narq4,narq6,narq9,narq10,narq11,narq12,narq13,narq14,narq17))
dat.im <- subset(dat1, select = c(im1r,im2,im3r,im4r,im5r,im6,im7,im8))
dat.pref_index <- subset(dat1, select = c(Zdescriptions,Zlbdq))
round(psych::alpha(dat.narq_adm, na.rm=T)$total[1], 2)
round(psych::alpha(dat.narq_riv, na.rm=T)$total[1], 2)
round(psych::alpha(dat.im, na.rm=T)$total[1], 2)
round(psych::alpha(dat.pref_index, na.rm=T)$total[1], 2)
rm(dat.narq_adm, dat.narq_riv, dat.im, dat.pref_index)


#Table 1: Descriptive Statistics and Zero-Order Correlations in Study 1
dat_table1 <- data.frame(matrix(nrow=4, ncol=3))
dat_table1[,1] <- c("ADM","RIV","IM","Leadership Preferences")
dat_table1[,2] <- c(mean(dat1$narq_adm, na.rm=T),mean(dat1$narq_riv, na.rm=T),mean(dat1$IM, na.rm=T),mean(dat1$pref_index, na.rm=T))
dat_table1[,3] <- c(sd(dat1$narq_adm, na.rm=T),sd(dat1$narq_riv, na.rm=T),sd(dat1$IM, na.rm=T),sd(dat1$pref_index, na.rm=T))
dat1_1 <- subset(dat1, select = c(narq_adm,narq_riv,IM,pref_index))
cor_dat1_1 <- corr.test(dat1_1)
dat_cor <- as.data.frame(matrix(unlist(cor_dat1_1), nrow=157, ncol = 1, byrow=F), stringsAsFactors=FALSE)
dat_cor_r <- as.data.frame(matrix(dat_cor[1:16,], nrow=4))
dat_cor_p <- as.data.frame(matrix(dat_cor[34:49,], nrow=4))

dat_table1[,2:3] <- sapply(dat_table1[,2:3], as.numeric)
format_msd <- function(x){
  formatC(x, format = "f", digits = 2)
}
dat_table1[,2:3] <- dat_table1[,2:3] %>% mutate_if(is.numeric, format_msd)

dat_cor_r[,1:4] <- sapply(dat_cor_r[,1:4], as.numeric)
format_r <- function(x){
  ifelse(abs(x) < 0.005,  f_num(unlist(x), digits=3), f_num(unlist(x), digits=2))
}
dat_cor_r[,1:4] <- dat_cor_r[,1:4] %>% mutate_if(is.numeric, format_r)

dat_cor_p[,1:4] <- sapply(dat_cor_p[,1:4], as.numeric)
format_p <- function(x){
  ifelse(abs(x) < 0.05 & abs(x) >= 0.01,  "*", ifelse(abs(x) < 0.01 & abs(x) >= 0.001,  "**", ifelse(abs(x) < 0.001,  "***", "")))
}
dat_cor_p[,1:4] <- dat_cor_p[,1:4] %>% mutate_if(is.numeric, format_p)

dat_table1$r1 <- paste0(dat_cor_r[,1],dat_cor_p[,1])
dat_table1$r2 <- paste0(dat_cor_r[,2],dat_cor_p[,2])
dat_table1$r3 <- paste0(dat_cor_r[,3],dat_cor_p[,3])

dat_table1[1,4:6] <- NA
dat_table1[2,5:6] <- NA
dat_table1[3,6] <- NA

dat_table1$blank1 <- NA
dat_table1 <- dat_table1[,c(1:3,7,4:6)]

col_keys <- c("X1","X2","X3","blank1","r1","r2","r3")
head1 <- c("Variable","M","SD","","1","2","3")
head <- data.frame(col_keys,head1, stringsAsFactors = FALSE)
rm(col_keys,head1)

tbl <- flextable(dat_table1)
tbl <- set_header_df(tbl, mapping=head, key="col_keys")
tbl <- hline_top(tbl, j=1:7, border=fp_border(width=2), part="header")
tbl <- hline(tbl, i=1, j=1:7, border=fp_border(width=1.2), part="header")
tbl <- flextable::font(tbl, fontname="Times", part="all")
tbl <- fontsize(tbl, size=12, part="all")
tbl <- align(tbl, align="center", part="all")
tbl <- align(tbl, j = c("X1"), align="left", part="body")
tbl <- italic(tbl, i=1, j=c("X2","X3"), part="header")
tbl <- width(tbl, j =~ X1, width=1.8)
tbl <- width(tbl, j =~ X2 + X3, width=.55)
tbl <- width(tbl, j =~ blank1, width=.1)
tbl

setwd(files_wd)
doc <- read_docx()
doc <- body_add_flextable(doc, value = tbl)
print(doc, target = "Table_1.docx")

rm(cor_dat1_1, dat_cor, dat_cor_p, dat_cor_r, dat_table1, dat1_1, doc, head, tbl, format_msd, format_p, format_r)


#Note of Table 1
round(mean(dat1$descriptions),2)
round(sd(dat1$descriptions),2)
round(mean(dat1$lbdq),2)
round(sd(dat1$lbdq),2)
t.test(dat1$descriptions, mu = 3.5, alternative = "two.sided")
cohens_d(dat1$descriptions, mu=3.5)
t.test(dat1$lbdq, mu = 3.5, alternative = "two.sided")
cohens_d(dat1$lbdq, mu=3.5)


#Table 2: Unique Relations of Narcissistic Admiration and Rivalry to Leadership Preferences in Study 1

##regression models
model1 <- lm(Zpref_index ~ Znarq_adm + Znarq_riv, data=dat1)
dat1_1 <- rbind(data.frame(summary(model1)$coef,confint(model1)))
model1_rsquared_p <- pf(summary(model1)$fstatistic[1],summary(model1)$fstatistic[2],summary(model1)$fstatistic[3],lower.tail=FALSE)
model1_rsquared_p <- ifelse(model1_rsquared_p < 0.001,  "< .001", paste0("= ", substr(round(model1_rsquared_p, 3),2,5)))
model1_rsquared <- paste0("R2 = ",substr(round(summary(model1)$r.squared,2),2,4),", F(",summary(model1)$fstatistic[2],", ",summary(model1)$fstatistic[3],") = ",round(summary(model1)$fstatistic[1],2),", p ",model1_rsquared_p)

model2 <- lm(Zpref_index ~ Znarq_adm + Znarq_riv + ZIM, data=dat1)
dat1_2 <- rbind(data.frame(summary(model2)$coef,confint(model2)))
model2_rsquared_p <- pf(summary(model2)$fstatistic[1],summary(model2)$fstatistic[2],summary(model2)$fstatistic[3],lower.tail=FALSE)
model2_rsquared_p <- ifelse(model2_rsquared_p < 0.001,  "< .001", paste0("= ", substr(round(model2_rsquared_p, 3),2,5)))
model2_rsquared <- paste0("R2 = ",substr(round(summary(model2)$r.squared,2),2,4),", F(",summary(model2)$fstatistic[2],", ",summary(model2)$fstatistic[3],") = ",round(summary(model2)$fstatistic[1],2),", p ",model2_rsquared_p)


##effect size of unique relation between narcissistic rivalry and leadership preferences
model1_no_riv <- lm(Zpref_index ~ Znarq_adm, data=dat1)
cohens_f_squared_riv <- cohens_f_squared(model1, model2 = model1_no_riv)
fstatistic_riv <- anova(model1,model1_no_riv)
fstatistic_riv_p <- ifelse(fstatistic_riv[2,6] < 0.001,  "< .001", paste0("= ", substr(round(fstatistic_riv[2,6], 3),2,5)))
paste0("f2 = ",substr(round(cohens_f_squared_riv[1,1],2),2,4),", ΔR2 = ",substr(round(cohens_f_squared_riv[1,5],2),2,4),", F(1, 313) = ",round(fstatistic_riv[2,5],2),", p ",fstatistic_riv_p)


##create table
dat1_1 <- dat1_1[!(row.names(dat1_1) %in% "(Intercept)"), ]
dat1_2 <- dat1_2[!(row.names(dat1_2) %in% "(Intercept)"), ]
dat1_1[nrow(dat1_1)+1,] <- NA
dat1_1$Predictor <- NA
dat1_1 <- dat1_1[,c(7,1,5:6,3:4)]
dat1_1$Predictor <- c("ADM","RIV","IM")
dat1_2$blank <- NA
dat_table2 <- cbind(dat1_1,dat1_2[,c(7,1,5:6,3:4)])

colnames(dat_table2) <- c("Predictor","beta1","ci_lower1","ci_upper1","t1","p1","blank","beta2","ci_lower2","ci_upper2","t2","p2")

dat_table2[,2:12] <- sapply(dat_table2[,2:12], as.numeric)

format_beta <- function(x){
  ifelse(abs(x) < 0.005,  f_num(unlist(x), digits=3), f_num(unlist(x), digits=2))
}
dat_table2[,c(2:4,8:10)] <- dat_table2[,c(2:4,8:10)] %>% mutate_if(is.numeric, format_beta)

format_t <- function(x){
  formatC(x, format = "f", digits = 2)
}
dat_table2[,c(5,11)] <- dat_table2[,c(5,11)] %>% mutate_if(is.numeric, format_t)

format_p <- function(x){
  ifelse(abs(x) < 0.001,  "< .001", f_num(unlist(x), digits=3))
}
dat_table2[,c(6,12)] <- dat_table2[,c(6,12)] %>% mutate_if(is.numeric, format_p)

dat_table2$ci1 <- paste0("[", dat_table2$ci_lower1, ", ", dat_table2$ci_upper1, "]")
dat_table2$ci2 <- paste0("[", dat_table2$ci_lower2, ", ", dat_table2$ci_upper2, "]")

dat_table2 <- dat_table2[,c("Predictor","beta1","ci1","t1","p1","blank","beta2","ci2","t2","p2")]
dat_table2[3,c(3:4)] <- NA

dat_table2[4,2] <- model1_rsquared
dat_table2[4,7] <- model2_rsquared

col_keys <- c("Predictor","beta1","ci1","t1","p1","blank","beta2","ci2","t2","p2")
head1 <- c("","Model 1","Model 1","Model 1","Model 1","","Model 2","Model 2","Model 2","Model 2")
head2 <- c("Predictor","\u03B2","95% CI","t","p","","\u03B2","95% CI","t","p")
head <- data.frame(col_keys,head1,head2, stringsAsFactors = FALSE)
rm(col_keys,head1,head2)

tbl <- flextable(dat_table2)
tbl <- set_header_df(tbl, mapping=head, key="col_keys")
tbl <- hline_top(tbl, j=1:10, border=fp_border(width=2), part="header")
tbl <- merge_at(tbl, i=1, j=2:5, part="header")
tbl <- merge_at(tbl, i=1, j=7:10, part="header")
tbl <- hline(tbl, i=1, j=c(2:5,7:10), border=fp_border(width=1.2), part="header")
tbl <- hline(tbl, i=2, j=1:10, border=fp_border(width=1.2), part="header")
tbl <- merge_at(tbl, i=4, j=2:5, part="body")
tbl <- merge_at(tbl, i=4, j=7:10, part="body")
tbl <- flextable::font(tbl, fontname="Times", part="all")
tbl <- fontsize(tbl, size=12, part="all")
tbl <- fontsize(tbl, i=4, size=11, part="body")
tbl <- align(tbl, align="center", part="all")
tbl <- align(tbl, j = c("Predictor"), align="left", part="body")
tbl <- italic(tbl, i=2, j=c("t1","p1","t2","p2"), part="header")
tbl <- width(tbl, j =~ Predictor, width=.9)
tbl <- width(tbl, j =~ beta1 + beta2, width=.45)
tbl <- width(tbl, j =~ ci1 + ci2, width=.95)
tbl <- width(tbl, j =~ t1 + t2 + p1 + p2, width=.6)
tbl <- width(tbl, j =~ blank, width=.3)
tbl <- colformat_lgl(tbl, j=2:6, na_str="")
tbl

setwd(files_wd)
doc <- read_docx()
doc <- body_add_flextable(doc, value = tbl)
print(doc, target = "Table_2.docx")

rm(model1, model2, model1_no_riv, fstatistic_riv, cohens_f_squared_riv, dat_table2, dat1_1, dat1_2, doc, head, tbl, fstatistic_riv_p, model1_rsquared, model1_rsquared_p, model2_rsquared, model2_rsquared_p, format_beta, format_p, format_t)


#Note of Table 2

##semi-partial correlation between impression management and leadership preferences in Model 2
spcor_IM <- spcor.test(dat1$pref_index,dat1$IM,dat1[,c("narq_adm","narq_riv")])
paste0("sr(312) = ",round(spcor_IM[1],2)," p = ",substr(round(spcor_IM[2],3),2,5))

##semi-partial correlation between impression management and leadership preferences when narcissistic admiration was omitted from Model 2
spcor_IM_RIV <- spcor.test(dat1$pref_index,dat1$IM,dat1$narq_riv)
paste0("sr(313) = ",round(spcor_IM_RIV[1],2)," p = ",substr(round(spcor_IM_RIV[2],3),2,5))

##semi-partial correlation between impression management and leadership preferences when narcissistic rivalry was omitted from Model 2
spcor_IM_ADM <- spcor.test(dat1$pref_index,dat1$IM,dat1$narq_adm)
paste0("sr(313) = ",round(spcor_IM_ADM[1],2)," p = ",substr(round(spcor_IM_ADM[2],3),2,5))

rm(spcor_IM, spcor_IM_RIV, spcor_IM_ADM)


#Footnote 5
model1a <- lm(descriptions ~ Znarq_adm + Znarq_riv, data=dat1)
model1a_coef <- data.frame(summary(model1a)$coef)
Znarq_riv_reversed_PDDA_descriptions <- (3.49 - model1a_coef[1,1]) / (model1a_coef[3,1])
model1b <- lm(lbdq ~ Znarq_adm + Znarq_riv, data=dat1)
model1b_coef <- data.frame(summary(model1b)$coef)
Znarq_riv_reversed_PDDA_lbdq <- (3.49 - model1b_coef[1,1]) / (model1b_coef[3,1])
paste0("The PDDA was reversed (i.e., preference for autocratic over democratic leadership) for narcissistic rivalry scores ≥ M + ",round(Znarq_riv_reversed_PDDA_descriptions,2)," SD (leadership descriptions measure) and M + ",round(Znarq_riv_reversed_PDDA_lbdq,2)," SD (leadership behaviors measure).")


#clear environment
rm(list = ls())