Sex differences
In this section, several different ways for estimating sex differences in personality are demonstrated. We start with typical univariate differences in the broad Big Five domains, followed by multivariate differences in both domains and using personality items as a multivariate set. Special attention is given to regularized methods because sex differences and the distributions underlying these mean differences can be both examined from this approach.
Univariate differences in trait domain sum scores using the d_pooled_sd -function
d_pooled_sd -function is used to calculate Cohen’s d with pooled standard deviation across the male and female groups for each Big Five trait separately. The arguments include the analyzed Big Five trait (var), sex as the binary grouping variable (group.var), the values that represents these groups (group.values), and logical for whether statistical inferences is wanted (infer). All arguments except statistical inference are also included in the multivariate -function D_regularized presented further below.
d.N<-d_pooled_sd(data = US.dat,var = "N",
group.var = "gender.cat",group.values = c("Male","Female"),
infer=T)
d.E<-d_pooled_sd(data = US.dat,var = "E",
group.var = "gender.cat",group.values = c("Male","Female"),
infer=T)
d.O<-d_pooled_sd(data = US.dat,var = "O",
group.var = "gender.cat",group.values = c("Male","Female"),
infer=T)
d.A<-d_pooled_sd(data = US.dat,var = "A",
group.var = "gender.cat",group.values = c("Male","Female"),
infer=T)
d.C<-d_pooled_sd(data = US.dat,var = "C",
group.var = "gender.cat",group.values = c("Male","Female"),
infer=T)
# compile to a table
uni.d<-rbind(d.N,d.E,d.O,d.A,d.C)
rownames(uni.d)<-c("N","E","O","A","C")
round(uni.d,2)#> n.Male n.Female m.Male m.Female sd.Male sd.Female pooled.sd diff D
#> N 2945 5739 2.79 3.16 0.88 0.86 0.87 -0.37 -0.42
#> E 2945 5739 3.00 3.06 0.96 0.96 0.96 -0.07 -0.07
#> O 2945 5739 4.07 3.88 0.58 0.63 0.62 0.18 0.30
#> A 2945 5739 3.68 4.03 0.75 0.68 0.71 -0.35 -0.50
#> C 2945 5739 3.38 3.44 0.72 0.75 0.74 -0.07 -0.09
#> t_stat df p
#> N -18.45 5845.38 0
#> E -3.09 5919.79 0
#> O 13.44 6458.37 0
#> A -21.33 5455.81 0
#> C -3.95 6119.69 0Column D indicates Cohen’s d. Females score higher in Neuroticism (N; d = -0.42) and Agreeableness (A; d = -0.50), while Males score higher in Openness (O; d = 0.30). Differences in Extraversion (E) and Conscientiousness (C), albeit statistically significant (p < .001), seem negligible in magnitude (d < |.10|) .
The average absolute difference across the Big Five traits is therefore: d = 0.27. Importantly, this average absolute difference is not informative about the overall distance between males and females. One should use multivariate methods for examining questions regarding the overall distance in certain domain, like personality.
Other columns in the output are:
n.Male and n.Female = sample sizes
m.Male and m.Female = mean levels
sd.Male and sd.Female = standard deviations
pooled.sd = pooled standard deviation
diff = the raw (non-standardized) difference
t_stat = t statistic
df = Welch modified degrees of freedom
p = p-value
Multivariate difference using trait sum scores as variables
Although the univariate differences presented above estimate sex difference in each broad Big Five trait separately, these single Cohen’s d estimates do not provide a clear picture about the overall distance between males and females. Multivariate methods, using multiple variables and differences in these variables at the same, while accounting for their covariability (intercorrelations), should be used for understanding overall distances. These D values can be interpreted in the same way as Cohen’s d, i.e., standardized distances between males and females.
Mahalanobis’ D
Mahalanobis’ D is the standardized distance between the Male and Female group centroids in multivariate space. For calculation of Mahalanobis’ D, correlation matrix of the multivariate set and vectors of standardized mean differences are needed. The ordering of the variables should be same in these inputs.
# let's use the previously obtained univariate d-values
d_vector<-uni.d[,"D"]
# calculate intercorrelations
cor_mat<-cor(US.dat[,c("N","E","O","A","C")])
# test that the ordering of these variables matches
names(d_vector)==colnames(cor_mat)#> [1] TRUE TRUE TRUE TRUE TRUE# print the correlation matrix
round(cor_mat,2)#> N E O A C
#> N 1.00 -0.25 -0.10 -0.10 -0.26
#> E -0.25 1.00 0.15 0.31 0.10
#> O -0.10 0.15 1.00 0.10 0.05
#> A -0.10 0.31 0.10 1.00 0.17
#> C -0.26 0.10 0.05 0.17 1.00# Calculate D
D_maha<-sqrt(t(d_vector) %*% solve(cor_mat) %*% d_vector)
D_maha#> [,1]
#> [1,] 0.7681502Mahalanobis’ D is, of course, larger than any of the univariate d’s for single variables, here the standardized distance between the group centroids of males and females in the multivariate space is D = 0.77
Regularized D with D_regularized -function
Regularized variant of D is calculated using a predictive approach whereby the binary sex-variable is predicted by the multivariate variable set with regularized logistic regression. The predicted values for each individual are indicative of logits of probabilities of being a male or female.
These femininity-masculinity scores (FM) can be used in similar ways to any other variable to calculate standardized distances between mean FM of Males and mean FM of females with pooled standard deviations. Beyond estimation of mean differences, there are many other interesting features that can be informative of the multivariate/FM distributions.
In the D_regularized -function, separate partitions of the original data are used for regulation (i.e., training) and prediction (i.e., testing). Regularization means that the beta-coefficients for variables in the multivariate set are penalized through a specified function (here, a default elastic net penalty with 10-fold cross-validation is used). After this, the obtained coefficients weights are used in the prediction part of the data as measures of FM for each individual. Argument out (logical, meaning are predictions made “out-of-bag”) and size (size of the training partition of the data, equal for males and females) are used for informing the function that a partition is wanted. The multivariate set is supplied with mv.vars -argument.
By requesting, it is possible to calculate pcc (probability of correctly classified), auc (area under the receiver operating characteristics), and also frequencies indicating the proportions of the sample that fall within certain cutoffs of predicted probabilities from 0 to 1 for each group (pred.prob).
# set seed number for reproducibility
set.seed(37127)
# decide the size of the regularization (training data fold)
# here, half of the smaller of the male and female groups is used
size=round(min(table(US.dat$gender.cat))/2)
# run D_regularized
D_out<-
D_regularized(data=US.dat,
mv.vars=c("N","E","O","A","C"),
group.var="gender.cat",
group.values = c("Male","Female"),
out = T,
size=size,
pcc = T,
auc = T,
pred.prob = T)
# print the output
round(D_out$D,2)#> n.Male n.Female m.Male m.Female sd.Male sd.Female pooled.sd diff D
#> [1,] 1473 4267 0.34 -0.3 0.86 0.79 0.81 0.64 0.79
#> pcc.Male pcc.Female pcc.total auc
#> [1,] 0.63 0.67 0.66 0.71From the output, we see that the sample sizes are smaller as compared to univariate estimates, because the sample sizes are only indicative of the samples in prediction/testing dataset.
Standardized multivariate difference between males and females was estimated at D = 0.79. Very similar in magnitude to multivariate sex difference indexed by Mahalanobis’ D.
There are additional features in this output that might be of interest. Similarly to univariate d estimates, we get descriptive statistics for both groups on FM-score metric, as well as an estimate of pooled estimate of standard deviation that is also used for calculation of D. pcc indicates probability of correctly classified individuals, and auc is the area under the receiver operating characteristics. These are both commonly used indices of predictive performance. In this exemplary case, the performance of broad Big Five domains in predicting sex. Detailing predicted probabilities of being a male across certain cutoffs (defaults were used in this example) for both groups can be obtained from P.table part of the output:
round(100*D_out$P.table,1)#>
#> [0,0.2) [0.2,0.4) [0.4,0.6) [0.6,0.8) [0.8,1]
#> Female 7.2 39.0 37.0 14.1 2.7
#> Male 1.2 18.1 37.1 32.6 11.0From this table we see that majority of males are nevertheless not very clearly classified as males based on the Big Five personality traits (11% in .80-1 range), and same is true for classifying females as females (7% in 0-0.20 range).
For examining the regularized beta weights obtained at the training phase, call coefficients on the cv.mod object of the output. Argument s is the value of the penalty parameter lambda by which the predictions are calculate. The default “lambda.min” uses the values where the mean cross-validated error is minimized, and “lambda.1se” uses the largest value of lambda such that error is within 1 standard error of the minimum. To obtain predicted values for a certain s, it must be supplied already to D_regularized -function, but the weights given by different approaches can be examined irregardless of the chosen value of s. It is also possible to plot the cross-validation curve to visualize how “lambda.min” and “lambda.1se” were obtained. For more information on regularized regression conducted with the glmnet package, please refer to https://glmnet.stanford.edu/articles/glmnet.html.
coefficients(D_out$cv.mod,
s = "lambda.min")#> 6 x 1 sparse Matrix of class "dgCMatrix"
#> s1
#> (Intercept) 3.37297381
#> N -0.54865910
#> E -0.02070834
#> O 0.64983681
#> A -0.88854735
#> C -0.23751055coefficients(D_out$cv.mod,
s = "lambda.1se")#> 6 x 1 sparse Matrix of class "dgCMatrix"
#> s1
#> (Intercept) 2.3670214
#> N -0.3900418
#> E .
#> O 0.4626928
#> A -0.6930196
#> C -0.1053781plot(D_out$cv.mod)
Overlap
The same pattern can also be observed from visually inspecting the distributions of the FM-scores, and assessing the overlap between the distributions of FM-scores of males and females.
# Obtain the D-estimate
D<-unname(D_out$D[,"D"])
# Overlap relative to a single distribution (OVL)
OVL<-2*pnorm((-D/2))
# Proportion of overlap relative to the joint distribution
OVL2<-OVL/(2-OVL)
# non-parametric overlap with overlapping package
np.overlap<-
overlap(x = list(D_out$pred.dat[
D_out$pred.dat$group=="Male","pred"],
D_out$pred.dat[
D_out$pred.dat$group=="Female","pred"]),
plot=T)
# this corresponds to Proportion of overlap relative to the joint distribution (OVL2)
np.OVL2<-unname(np.overlap$OV)
# from which Proportion of overlap relative to a single distribution (OVL) is approximated at
np.OVL<-(2*np.OVL2)/(1+np.OVL2)
# comparison of different overlap-estimates
round(cbind(OVL,np.OVL,OVL2,np.OVL2),2)#> OVL np.OVL OVL2 np.OVL2
#> [1,] 0.69 0.7 0.53 0.54In this example, the non-parametric overlap (np) estimates to these logistic distributions are almost identical to the parametric overlap estimates.
Bootstrap confidence intervals for Regularized D
To obtain confidence intervals for D-estimates in the regularized techniques, a bootstrap approach can be used whereby D is estimated with multiple redrawn bootstrap samples from the original sample and using the distribution of these D values for constructing a percentile confidence interval around the estimate. In this example, we take only 10 redraws with replacement, but in reality one should take at least 1000 or 5000.
# number of redraws
reps=10
# placeholder for D-estimates
D_out_boot<-rep(NA,reps)
# repeat the procedure to obtain a distribution
t1<-Sys.time()
for (i in 1:reps){
# draw replaced samples with the same size as the original data
boot.dat<-sample_n(US.dat,size=nrow(US.dat),replace=T)
# run D_regularized for each sample
D_out_boot[i]<-
D_regularized(data=boot.dat,
mv.vars=c("N","E","O","A","C"),
group.var="gender.cat",
group.values = c("Male","Female"),
out = T,
size=size)$D[,"D"]
}
t2<-Sys.time()
t2-t1#> Time difference of 6.255088 secs# print 95% CI percentile interval
round(quantile(D_out_boot,c(.025,.975)),2)#> 2.5% 97.5%
#> 0.77 0.83Thus, it can be reported that the multivariate sex difference D = 0.79, 95% CI [0.77, 0.83].
Multivariate difference using item scores as variables
Recently there have been multiple studies in the field of personality psychology highlighting the importance of narrow personality traits, also in terms of their ability to predict criteria over and above the more commonly used broad trait domains, such as Big Five (e.g., Seeboth & Mõttus, 2018). These narrow personality traits are usually either facets or nuances that are commonly modeled as being part of a certain Big Five trait, although their unique predictive utility suggests that they have trait-like properties (stability, cross-informant agreement etc.) that go beyond the broad Big Five domains (e.g., Mõttus et al., 2019). That is, the narrower traits can be considered unique traits as well. It is therefore interesting to see whether the item scores as the variable set in multivariate estimation would also better predictions of sex, and different estimates of sex differences (see Ilmarinen et al., 2022).
Because the number of the variables is increased (from five to fifty in this dataset), the use of regularization methods and different data partitions for training and testing are strongly recommended to avoid overfitting the data (for illustrative purposes, Mahalanobis’ D is also be calculated below).
Mahalanobis’ D
# place-holder vector for univariate Cohen d's for items
d_vector_items<-rep(NA,length(per.items))
# loop through the items and obtain d for each item
for (i in 1:length(per.items)){
d_vector_items[i]<-
d_pooled_sd(data=US.dat,
var=per.items[i],
group.var="gender.cat",
group.values=c("Male","Female"))[,"D"]
}
cor_mat_items<-cor(US.dat[,per.items])
# Calculate D
D_maha_items<-
sqrt(t(d_vector_items) %*% solve(cor_mat_items) %*% d_vector_items)
D_maha_items#> [,1]
#> [1,] 0.9212846Item-based Mahalanobis’ D is somewhat larger that any of those based on Big Five domains. The standardized distance between the group centroids of males and females in a space comprised of personality items is D = 0.92
Regularized D with D_regularized -function
It is straightforward to use D_regularized -function even with large variable sets (especially when the variable names are save in a vector first). Using separate fold for regularization and prediction for items. Domain-based results are printed for comparison.
D_items_out<-
D_regularized(data=US.dat,
mv.vars=per.items,
group.var="gender.cat",
group.values = c("Male","Female"),
out = T,
size=size,pcc = T,auc = T,pred.prob = T)
# Item-based D
round(D_items_out$D,2)#> n.Male n.Female m.Male m.Female sd.Male sd.Female pooled.sd diff D
#> [1,] 1473 4267 0.46 -0.4 0.9 0.86 0.87 0.86 0.99
#> pcc.Male pcc.Female pcc.total auc
#> [1,] 0.7 0.69 0.69 0.76# Big Five domain-based D
round(D_out$D,2)#> n.Male n.Female m.Male m.Female sd.Male sd.Female pooled.sd diff D
#> [1,] 1473 4267 0.34 -0.3 0.86 0.79 0.81 0.64 0.79
#> pcc.Male pcc.Female pcc.total auc
#> [1,] 0.63 0.67 0.66 0.71With items as the multivariate set, D = 0.99 is slightly larger than with Big Five trait domains, D = 0.79.
Differences in predictive performance can also be illustrated from comparing the predicted probabilities.
round(100*D_items_out$P.table,1)#>
#> [0,0.2) [0.2,0.4) [0.4,0.6) [0.6,0.8) [0.8,1]
#> Female 11.6 38.5 33.4 14.3 2.2
#> Male 1.6 15.1 32.9 35.4 15.0round(100*D_out$P.table,1)#>
#> [0,0.2) [0.2,0.4) [0.4,0.6) [0.6,0.8) [0.8,1]
#> Female 7.2 39.0 37.0 14.1 2.7
#> Male 1.2 18.1 37.1 32.6 11.0We see that there are more people in the extremes, but based on the pcc and auc -values printed above, the overall predictive performance is not drastically improved. A large proportion (more than 30%) of the data is still classified in the 40%-60% range, meaning their personality profiles are rather indifferent in regards to sex.
Overlap
Overlaps give similar information but in another format:
# Obtain the D-estimate
D_items<-unname(D_items_out$D[,"D"])
# Overlap relative to a single distribution (OVL)
OVL_items<-2*pnorm((-D_items/2))
# Proportion of overlap relative to the joint distribution
OVL2_items<-OVL_items/(2-OVL_items)
# non-parametric overlap
np.overlap_items<-
overlap(x = list(D_items_out$pred.dat[
D_items_out$pred.dat$group=="Male","pred"],
D_items_out$pred.dat[
D_items_out$pred.dat$group=="Female","pred"]),
plot=T)
# this corresponds to Proportion of overlap relative to the joint distribution (OVL2)
np.OVL2_items<-unname(np.overlap_items$OV)
# from which Proportion of overlap relative to a single distribution (OVL) is approximated at
np.OVL_items<-(2*np.OVL2_items)/(1+np.OVL2_items)
# compare the different overlap-estimates
round(cbind(OVL_items,np.OVL_items,
OVL2_items,np.OVL2_items),2)#> OVL_items np.OVL_items OVL2_items np.OVL2_items
#> [1,] 0.62 0.63 0.45 0.46# print Big Five domain based overlaps for comparison
round(cbind(OVL,np.OVL,OVL2,np.OVL2),2)#> OVL np.OVL OVL2 np.OVL2
#> [1,] 0.69 0.7 0.53 0.54There overlap between FM-distributions of males and females is somewhat smaller when individual items are used as the multivariate set rather than broad Big Five domains.