library(haven)
library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.0.5
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
AGagingData_2 <- read_sav("C:/Users/kimng/Desktop/Aging, Emotional Recognition, and Amygdala/AGagingData-2.sav")

# ++++++++++++++++++++++++++++
# flattenCorrMatrix
# ++++++++++++++++++++++++++++
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    cor  =(cormat)[ut],
    p = pmat[ut]
    )
}

3.1. P-values of pairwise correlations between Age and ER-40 scores and scores for each emotion

With Whole sample

df <- AGagingData_2[, c('Age','FEST', 'Neutral', 'Happy', 'Sad', 'Angry', 'Fearful')]
res<-rcorr(as.matrix(df)) #correlation matrix
flat_corr <- flattenCorrMatrix(res$r, res$P)
age_corr <- flat_corr[flat_corr$row == "Age" ,] #get p values for correlation with age only
age_corr$adj_p <- p.adjust(age_corr$p, 'fdr') # FDR adjusted p-values
age_corr #Correlation table of correlations between age and ER-40 performance
##    row  column         cor            p        adj_p
## 1  Age    FEST -0.33629608 0.0001180941 0.0007085649
## 2  Age Neutral -0.26180455 0.0030635179 0.0061270358
## 4  Age   Happy -0.28745835 0.0010995846 0.0032987539
## 7  Age     Sad -0.02753025 0.7596012310 0.7596012310
## 11 Age   Angry -0.24318318 0.0060725845 0.0091088768
## 16 Age Fearful -0.21147705 0.0174485605 0.0209382726

Analyses by gender

#Male
df_m <- AGagingData_2[AGagingData_2$Sex ==1, c('Age','FEST', 'Neutral', 'Happy', 'Sad', 'Angry', 'Fearful')]
res_m <-rcorr(as.matrix(df_m)) #correlation matrix
flat_corr_m <- flattenCorrMatrix(res_m$r, res_m$P)
age_corr_m <- flat_corr_m[flat_corr_m$row == "Age" ,] #get p values for correlation with age only
age_corr_m <- cbind(gender = 'Male', age_corr_m)

#Fenale
df_f <- AGagingData_2[AGagingData_2$Sex ==2, c('Age','FEST', 'Neutral', 'Happy', 'Sad', 'Angry', 'Fearful')]
res_f <-rcorr(as.matrix(df_f)) #correlation matrix
flat_corr_f <- flattenCorrMatrix(res_f$r, res_f$P)
age_corr_f <- flat_corr_f[flat_corr_f$row == "Age" ,] #get p values for correlation with age only
age_corr_f <- cbind(gender = 'Female', age_corr_f)


age_corr_all <- rbind(age_corr_m, age_corr_f) # combining correlation tables for males and females 
age_corr_all$adj_p <- p.adjust(age_corr_all$p, 'fdr') # FDR adjusted p-values
age_corr_all
##     gender row  column         cor            p       adj_p
## 1     Male Age    FEST -0.28304344 0.0313259516 0.075182284
## 2     Male Age Neutral -0.09128358 0.4955637507 0.594676501
## 4     Male Age   Happy -0.23644343 0.0739484225 0.126768724
## 7     Male Age     Sad -0.08001385 0.5504668369 0.600509277
## 11    Male Age   Angry -0.10571824 0.4296349271 0.572846569
## 16    Male Age Fearful -0.27116036 0.0395047959 0.079009592
## 12  Female Age    FEST -0.38116242 0.0013419079 0.007249295
## 21  Female Age Neutral -0.40710319 0.0005703222 0.006843867
## 41  Female Age   Happy -0.33593255 0.0050991784 0.015297535
## 71  Female Age     Sad  0.01770623 0.8860426070 0.886042607
## 111 Female Age   Angry -0.37152913 0.0018123237 0.007249295
## 161 Female Age Fearful -0.15654170 0.2023731757 0.303559764

3.2 P-values of pairwise correlations between age and miscategorization errors

Exploratory analyses with other types of errors in results section 3.2

Males

df_m <- AGagingData_2[AGagingData_2$Sex ==1 
                      , c('Age','Neutral_NH', 'Neutral_NS', 'Neutral_NA', 'Neutral_NF',
                          'Sad_SN', 'Sad_SH', 'Sad_SA', 'Sad_SF',
                          'Angry_AN', 'Angry_AH', 'Angry_AS', 'Angry_AF'
                          )]
res_m <-rcorr(as.matrix(df_m)) #correlation matrix
flat_corr_m <- flattenCorrMatrix(res_m$r, res_m$P)
age_corr_m <- flat_corr_m[flat_corr_m$row == "Age" ,]
age_corr_m$adj_p <- p.adjust(age_corr_m$p, 'fdr') # FDR adjusted p-values
age_corr_m
##    row     column         cor          p     adj_p
## 1  Age Neutral_NH -0.26715648 0.04262932 0.2712482
## 2  Age Neutral_NS  0.08687653 0.51668952 0.5166895
## 4  Age Neutral_NA  0.18619875 0.16167788 0.3233558
## 7  Age Neutral_NF  0.25919106 0.04945044 0.2712482
## 11 Age     Sad_SN  0.16940846 0.20362092 0.3292894
## 16 Age     Sad_SH -0.19925614 0.13373422 0.3233558
## 22 Age     Sad_SA -0.11392896 0.39446570 0.4303262
## 29 Age     Sad_SF  0.14622037 0.27341183 0.3645491
## 37 Age   Angry_AN  0.18661170 0.16073316 0.3233558
## 46 Age   Angry_AH -0.24149305 0.06781206 0.2712482
## 56 Age   Angry_AS  0.16368986 0.21952625 0.3292894
## 67 Age   Angry_AF -0.11992391 0.36989315 0.4303262

Females

df_f <- AGagingData_2[AGagingData_2$Sex ==2 
                      , c('Age', 'Sad_SN', 'Sad_SH', 'Sad_SA', 'Sad_SF',
                          'Fearful_FN', 'Fearful_FH', 'Fearful_FS', 'Fearful_FA'
                          )]

res_f<-rcorr(as.matrix(df_f)) #correlation matrix
flat_corr_f <- flattenCorrMatrix(res_f$r, res_f$P)
age_corr_f <- flat_corr_f[flat_corr_f$row == "Age" ,]
age_corr_f$adj_p <- p.adjust(age_corr_f$p, 'fdr') # FDR adjusted p-values
age_corr_f
##    row     column         cor          p     adj_p
## 1  Age     Sad_SN -0.01407284 0.90931574 0.9093157
## 2  Age     Sad_SH -0.07014698 0.56974750 0.7596633
## 4  Age     Sad_SA -0.03194286 0.79595156 0.9093157
## 7  Age     Sad_SF  0.09404261 0.44557644 0.7129223
## 11 Age Fearful_FN  0.10823069 0.37965562 0.7129223
## 16 Age Fearful_FH  0.09797557 0.42669293 0.7129223
## 22 Age Fearful_FS -0.11574746 0.34724502 0.7129223
## 29 Age Fearful_FA  0.27791008 0.02175484 0.1740387

3.3 Alexithymia traits, theory of mind, and emotional recognition

Relationship between age and TAS, RMET in males and females

df_m <- AGagingData_2[AGagingData_2$Sex ==1 , c('Age', 'TAS', 'RMET')] #get male data
df_f <- AGagingData_2[AGagingData_2$Sex ==2 , c('Age','TAS', 'RMET')] #get female data

res_m <-rcorr(as.matrix(df_m)) #correlation matrix
flat_corr_m <- flattenCorrMatrix(res_m$r, res_m$P)
age_corr_m <- flat_corr_m[flat_corr_m$row == "Age" ,] #get p values for correlation with age only
age_corr_m <- cbind(gender = 'Male', age_corr_m)

res_f <-rcorr(as.matrix(df_f)) #correlation matrix
flat_corr_f <- flattenCorrMatrix(res_f$r, res_f$P)
age_corr_f <- flat_corr_f[flat_corr_f$row == "Age" ,] #get p values for correlation with age only
age_corr_f <- cbind(gender = 'Female', age_corr_f)

age_corr_all <- rbind(age_corr_m, age_corr_f) # combining correlation tables for males and females 
age_corr_all$adj_p <- p.adjust(age_corr_all$p, 'fdr') # FDR adjusted p-values
age_corr_all
##   gender row column         cor         p     adj_p
## 1   Male Age    TAS -0.04331439 0.7468140 0.8800156
## 2   Male Age   RMET -0.21359456 0.1074161 0.4296644
## 3 Female Age    TAS -0.01864948 0.8800156 0.8800156
## 4 Female Age   RMET -0.11854503 0.3356357 0.6712714

Relationship between TAS, RMET and ER-40 performance in the entire sample

df <- AGagingData_2[, c('TAS', 'RMET', 'FEST', 'Neutral', 'Happy', 'Sad', 'Angry', 'Fearful')]
res<-rcorr(as.matrix(df)) #correlation matrix
flat_corr <- flattenCorrMatrix(res$r, res$P)
corr <- flat_corr[flat_corr$row == c("TAS")|flat_corr$row == c("RMET") ,] #get correlation with RMET and TAS only
corr <- corr[-1, ]
corr$adj_p <- p.adjust(corr$p, 'fdr')
corr
##     row  column          cor          p     adj_p
## 2   TAS    FEST -0.019957947 0.82445651 0.9634451
## 3  RMET    FEST  0.200257238 0.02455677 0.1591681
## 4   TAS Neutral  0.004124019 0.96344509 0.9634451
## 5  RMET Neutral  0.197646167 0.02652801 0.1591681
## 7   TAS   Happy  0.064831912 0.47076008 0.8070173
## 8  RMET   Happy  0.132635740 0.13872812 0.4161844
## 11  TAS     Sad -0.161930761 0.07005910 0.2802364
## 12 RMET     Sad  0.051425708 0.56740481 0.8511072
## 16  TAS   Angry -0.008091401 0.92834918 0.9634451
## 17 RMET   Angry  0.012565447 0.88893923 0.9634451
## 22  TAS Fearful  0.084423084 0.34727831 0.7327524
## 23 RMET Fearful  0.081146681 0.36637619 0.7327524

Following up with significant correlations in the entire sample

df_TAS <- AGagingData_2[, c('TAS', 'Sad_SN', 'Sad_SH', 'Sad_SA', 'Sad_SF')]
res_TAS<-rcorr(as.matrix(df_TAS)) #correlation matrix
flat_corr_TAS <- flattenCorrMatrix(res_TAS$r, res_TAS$P)
corr_TAS <- flat_corr_TAS[flat_corr_TAS$row == c("TAS"),] #get correlation TAS
df_RMET <- AGagingData_2[, c('RMET', 'Neutral_NH', 'Neutral_NS', 'Neutral_NA', 'Neutral_NF')]
res_RMET<-rcorr(as.matrix(df_RMET)) #correlation matrix
flat_corr_RMET <- flattenCorrMatrix(res_RMET$r, res_RMET$P)
corr_RMET <- flat_corr_RMET[flat_corr_RMET$row == c("RMET"),] #get correlation with RMET

corr_all <- rbind(corr_TAS, corr_RMET) # combining correlation tables 
corr_all$adj_p <- p.adjust(corr_all$p, 'fdr') # FDR adjusted p-values
corr_all
##     row     column          cor            p        adj_p
## 1   TAS     Sad_SN  0.198037263 2.622442e-02 0.0699317853
## 2   TAS     Sad_SH  0.161226776 7.130221e-02 0.1426044187
## 4   TAS     Sad_SA -0.001252541 9.888941e-01 0.9888941092
## 7   TAS     Sad_SF -0.008587947 9.239650e-01 0.9888941092
## 11 RMET Neutral_NH  0.019343879 8.297746e-01 0.9888941092
## 21 RMET Neutral_NS -0.209379076 1.862290e-02 0.0699317853
## 41 RMET Neutral_NA -0.038121654 6.717088e-01 0.9888941092
## 71 RMET Neutral_NF -0.370748365 1.924426e-05 0.0001539541

Following up with significant correlations between TAS and Sad recognition in males and females

df_m <- AGagingData_2[AGagingData_2$Sex ==1 , c('TAS', 'Sad_SN', 'Sad_SH')] #get male data
df_f <- AGagingData_2[AGagingData_2$Sex ==2 , c('TAS', 'Sad_SN', 'Sad_SH')] #get female data

res_m <-rcorr(as.matrix(df_m)) #correlation matrix
flat_corr_m <- flattenCorrMatrix(res_m$r, res_m$P)
age_corr_m <- flat_corr_m[flat_corr_m$row == "TAS" ,] #get p values for correlation with age only
age_corr_m <- cbind(gender = 'Male', age_corr_m)

res_f <-rcorr(as.matrix(df_f)) #correlation matrix
flat_corr_f <- flattenCorrMatrix(res_f$r, res_f$P)
age_corr_f <- flat_corr_f[flat_corr_f$row == "TAS" ,] #get p values for correlation with age only
age_corr_f <- cbind(gender = 'Female', age_corr_f)

age_corr_all <- rbind(age_corr_m, age_corr_f) # combining correlation tables for males and females 
age_corr_all$adj_p <- p.adjust(age_corr_all$p, 'fdr') # FDR adjusted p-values
age_corr_all
##   gender row column         cor           p      adj_p
## 1   Male TAS Sad_SN  0.29327388 0.025465616 0.05093123
## 2   Male TAS Sad_SH -0.13181624 0.323965445 0.43195393
## 3 Female TAS Sad_SN  0.07106661 0.564686617 0.56468662
## 4 Female TAS Sad_SH  0.34241921 0.004259662 0.01703865

Relationship between RMET and ER-40 in males and females

df_m <- AGagingData_2[AGagingData_2$Sex ==1 , c('RMET', 'FEST', 'Neutral', 'Happy', 'Sad', 'Angry', 'Fearful')] #get male data
df_f <- AGagingData_2[AGagingData_2$Sex ==2 , c('RMET', 'FEST', 'Neutral', 'Happy', 'Sad', 'Angry', 'Fearful')] #get female data

res_m <-rcorr(as.matrix(df_m)) #correlation matrix
flat_corr_m <- flattenCorrMatrix(res_m$r, res_m$P)
age_corr_m <- flat_corr_m[flat_corr_m$row == "RMET" ,] #get p values for correlation with age only
age_corr_m <- cbind(gender = 'Male', age_corr_m)

res_f <-rcorr(as.matrix(df_f)) #correlation matrix
flat_corr_f <- flattenCorrMatrix(res_f$r, res_f$P)
age_corr_f <- flat_corr_f[flat_corr_f$row == "RMET" ,] #get p values for correlation with age only
age_corr_f <- cbind(gender = 'Female', age_corr_f)

age_corr_all <- rbind(age_corr_m, age_corr_f) # combining correlation tables for males and females 
age_corr_all$adj_p <- p.adjust(age_corr_all$p, 'fdr') # FDR adjusted p-values
age_corr_all
##     gender  row  column          cor            p        adj_p
## 1     Male RMET    FEST  0.441629187 0.0005196980 0.0031181879
## 2     Male RMET Neutral  0.509857767 0.0000434957 0.0005219485
## 4     Male RMET   Happy  0.282219390 0.0318433736 0.1273734942
## 7     Male RMET     Sad  0.103781712 0.4381805472 0.9599468223
## 11    Male RMET   Angry  0.085777824 0.5220261120 0.9599468223
## 16    Male RMET Fearful  0.117970932 0.3777945105 0.9599468223
## 12  Female RMET    FEST  0.044275508 0.7199601167 0.9599468223
## 21  Female RMET Neutral  0.008223397 0.9469353728 0.9822380162
## 41  Female RMET   Happy  0.052920962 0.6682090612 0.9599468223
## 71  Female RMET     Sad -0.002750803 0.9822380162 0.9822380162
## 111 Female RMET   Angry -0.065875799 0.5935250513 0.9599468223
## 161 Female RMET Fearful  0.019884471 0.8721348139 0.9822380162

Following up with neutral and happy facial recognition errors and RMET in males

df <- AGagingData_2[AGagingData_2$Sex ==1 
                      , c('RMET', 'Neutral_NH', 'Neutral_NS', 'Neutral_NA', 'Neutral_NF',
                          'Happy_HN', 'Happy_HS', 'Happy_HA', 'Happy_HF')]
res<-rcorr(as.matrix(df)) #correlation matrix
flat_corr <- flattenCorrMatrix(res$r, res$P)
corr <- flat_corr[flat_corr$row == c("RMET") ,] #get correlation with RMET
corr$adj_p <- p.adjust(corr$p, 'fdr')
corr
##     row     column         cor            p        adj_p
## 1  RMET Neutral_NH  0.07054705 5.987277e-01 0.6842602165
## 2  RMET Neutral_NS -0.48998946 9.459908e-05 0.0003056008
## 4  RMET Neutral_NA -0.25524133 5.314925e-02 0.1062984960
## 7  RMET Neutral_NF -0.48489049 1.146003e-04 0.0003056008
## 11 RMET   Happy_HN -0.13474137 3.132484e-01 0.4176644850
## 16 RMET   Happy_HS  0.16774588 2.081585e-01 0.3330535554
## 22 RMET   Happy_HA -0.49501148 7.808189e-05 0.0003056008
## 29 RMET   Happy_HF -0.03600505 7.884498e-01 0.7884497969