Dear Ben,
I did some further work on the replication of one of the figures published by:
Del Giudice, M. (2023). The S-index: Summarizing patterns of sex differences at the distribution extremes. Personality and Individual Differences, 205, 112088. https://doi.org/10.1016/j.paid.2023.112088
The objective of the so-called S-index, proposed by Marco Del Giudice, is 'to compare groups that are about equally represented in the population, which makes it ideally suited to the analysis of [on] sex differences' but not limited to that.
You are in a better position than me to assess the mathematical and statistical appropriateness of what Marco proposes.
But, for my research, which does meet the relevant assumptions that he posits, it might be a useful addition to the set of tools for distributional analysis and comparison between groups, like females and males.
So, using dstat, I have worked on his example data and can now replicate one of the traits of Figure 1 (Openness).
I suppose that my code is something of a 'hack' that you might expect from a not too experienced coder like me.
Therefore, to replicate his figure:
using coefplot (somewhat improving it in the process) I go through the following:
*! https://www.statalist.org/forums/forum/general-stata-discussion/general/1429589-how-to-modify-orientation-of-y-title-text
*! inspired by splitvallabels.ado v1.1.2 Nicholas Winter/Ben Jann 14aug2008
*! Using your own code, but modified by Chen, 2019 March 5:
program define splitvarlabels2, rclass
version 8.2
syntax varlist(max=1) [if] [in], [ Length(int 2) local(string) nobreak]
local labname : variable label `varlist'
if "`break'"!="" {
local i 1
local maxi : word count `labname'
while `i'<=`maxi' {
if "`labname'"!="" {
local labpart : piece `i' `length' of `"`labname'"', nobreak
local newstring `"`newstring'`"`labpart'"' "'
local i=`i'+1
}
}
}
else {
local i 1
local maxi : length local labname
while `i'<=`maxi'/`length' {
if "`labname'"!="" {
local labpart : piece `i' `length' of `"`labname'"'
local newstring `"`newstring'`"`labpart'"' "'
local i=`i'+1
}
}
}
return local retitle `"`newstring'"'
if "`local'"!="" {
c_local `local' `"`newstring'"'
}
end
The paper and the dta file that I created using his R-code and data, I provide here in a zipped file:
PID 2023 Guidice - R dta.zip
Next, I use this code:
* FULL SCALE REPLICATION
use "ipip_domains_rep.dta", clear
////////////////////////////////////////////////////////////////////////////////////
* Run dstat Gaussian
local binary female
local var O
qui dstat density `var' , over(`binary') n(200) k(gau) ll(1) ul(5) nograph // if `constraint'
est sto gau
////////////////////////////////////////////////////////////////////////////////////
* Collect beta's and calculate their fraction
mat EB = e(b)'
mat EB0 = EB[1..200,1]
mat EB1 = EB[201..400,1]
mat EBT = EB0+EB1
* https://www.statalist.org/forums/forum/general-stata-discussion/general/1354470-divide-a-matrix
mata: st_matrix("EBf", st_matrix("EB0") :/ st_matrix("EBT"))
mata: st_matrix("EBm", st_matrix("EB1") :/ st_matrix("EBT"))
mat EBT = EBm+EBf
////////////////////////////////////////////////////////////////////////////////////
* Collect scale measurement points (at)
mat AT = e(at)'
mat AT0 = AT[1..200,1]
mat AT1 = AT[201..400,1]
* Only one is required
mat EBmf = AT0,EBm,EBf
matrix colnames EBmf = "patx" "patm" "patf"
svmat EBmf , names(pat)
tokenize patx_g patm_g patf_g
forvalues i = 1(1)3 {
ren pat`i' `1'
macro shift
}
* Data setup
gen patMg=patm_g*100
gen patFg=patf_g*100
label var patMg "%"
label var patFg "%"
* Get Y title-text to be rotated vertically using the program splitvarlabels2 (run before above)
splitvarlabels2 patFg, local(change)
local binary female
local var O
local vart "Openness by `binary' (Gaussian kernel function)"
qui sum `var'
local Llsd = r(mean)-(r(sd))
local L2sd = r(mean)-(2*r(sd))
local L3sd = r(mean)-(3*r(sd))
local L4sd = r(mean)-(4*r(sd))
local Ulsd = r(mean)+(r(sd))
local U2sd = r(mean)+(2*r(sd))
local U3sd = r(mean)+(3*r(sd))
local pM = r(mean)
local pL = r(mean)-(3*r(sd))
local pU = r(mean)+(3*r(sd))
coefplot (gau, keep(0:) ci(ci) ciopts(col(%30) lw(0)) axis(1) pstyle(p1)) ///
(gau, keep(1:) ci(ci) ciopts(col(maroon%30) lw(0)) axis(1) pstyle(p2)) ///
, at(at) recast(line) cirecast(rarea) graphreg(c(gs15) m(t+3 l-3 r-0)) ///
plotreg(m(l+0 r+1 t+15)) ytit("density", m(l+1 r+0)) name(d`var', replace) ///
ysc(range(0 1) noex) ylab(, labs(*1) angle(none) nogrid) ytit("Density", m(r-1) s(*1.2)) ///
xsc(range(1 5) noex) xmt(1.1(.1)4.9) xlabel(1(.5)5, nogrid) /// legend(off) ///
xtit("`vart'", m(t+0 b-2) s(*1.2)) legend(order(4 "females" 2 "males") ring(0) pos(11) region(lc(none)) symx(*.8)) ///
xline(`Llsd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`L2sd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`L3sd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`L4sd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`Ulsd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`U2sd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`U3sd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`pL' , lc(gs10) lw(*1) lp(shortdash)) /// S-index lower bound (selected)
xline(`pM' , lc(gs10) lw(*1) lp(dash)) /// S-index population mean
xline(`pU' , lc(gs10) lw(*1) lp(shortdash)) /// S-index upper bound (selected)
text(1.07 `pM' "M") ///
text(1.07 `L4sd' "-4 SD") text(1.07 `L3sd' "-3 SD") text(1.07 `L2sd' "-2 SD") text(1.07 `Llsd' "-1 SD") ///
text(1.07 `U3sd' "3 SD") text(1.07 `U2sd' "2 SD") text(1.07 `Ulsd' "1 SD")
addplot: (line patFg patx_g if patx_g>=`L3sd', yaxis(2) xaxis(2) lc(stblue)) (line patMg patx_g if patx_g>=`L3sd', yaxis(2) xaxis(2) lc(stred)), ysc(range(-100 70) alt axis(2) noex) xsc(range(1 5) alt axis(2) noex off) ytit(`change', orient(horizontal) m(l-6 r-2) s(*1.2) axis(2)) xtit("", axis(2)) legend(order(4 "females" 2 "males") ring(0) pos(11) region(lc(none)) symx(*.8))
which results in:
My opinion is that this plot is more clear, because of the visualization of the two distribution curves separate from the two 'patterns of proportions' with their scales. Also more informative is that dstat allows for the visualization of confidence intervals which allows for a further analysis of the distributional differences between the two groups.
Also note that Marco applies the Gaussian distribution function which, I suppose, follows from his assumption that the scale data under inspection more or less originates from an 'idealized' normal distribution. That might be the case, or not, and I also (tried to) advance his result but using the Epanechnikov distribution function from dstat, like:
* FULL SCALE Advancement
use "ipip_domains_rep.dta", clear
////////////////////////////////////////////////////////////////////////////////////
* Run dstat Epanechnikov
local var O
local vart Openness
local binary female
*return list
qui dstat density `var' , over(`binary') n(200) k(e) ll(1) ul(5) nograph
est sto epa
////////////////////////////////////////////////////////////////////////////////////
* Collect beta's and calculate their fraction
mat EB = e(b)'
mat EB0 = EB[1..200,1]
mat EB1 = EB[201..400,1]
mat EBT = EB0+EB1
* https://www.statalist.org/forums/forum/general-stata-discussion/general/1354470-divide-a-matrix
mata: st_matrix("EBf", st_matrix("EB0") :/ st_matrix("EBT"))
mata: st_matrix("EBm", st_matrix("EB1") :/ st_matrix("EBT"))
mat EBT = EBm+EBf
////////////////////////////////////////////////////////////////////////////////////
* Collect scale measurement points (at)
mat AT = e(at)'
mat AT0 = AT[1..200,1]
mat AT1 = AT[201..400,1]
* Only one is required
mat EBmf = AT0,EBm,EBf
matrix colnames EBmf = "patx" "patm" "patf"
svmat EBmf , names(pat)
tokenize patx_e patm_e patf_e
forvalues i = 1(1)3 {
ren pat`i' `1'
macro shift
}
////////////////////////////////////////////////////////////////////////////////////
* Data setup for coefplot (run before below)
gen patMe=patm_e*100
gen patFe=patf_e*100
label var patMe "%"
label var patFe "%"
* Get Y title-text to be rotated vertically using the program splitvarlabels2 (run before above)
splitvarlabels2 patFe, local(change)
local binary female
local var O
local vart "Openness by `binary' (Epanechnikov kernel function)"
qui sum `var'
local Llsd = r(mean)-(r(sd))
local L2sd = r(mean)-(2*r(sd))
local L3sd = r(mean)-(3*r(sd))
local L4sd = r(mean)-(4*r(sd))
local Ulsd = r(mean)+(r(sd))
local U2sd = r(mean)+(2*r(sd))
local U3sd = r(mean)+(3*r(sd))
local pM = r(mean)
local pL = r(mean)-(3*r(sd))
local pU = r(mean)+(3*r(sd))
coefplot (epa, keep(0:) ci(ci) ciopts(col(%30) lw(0)) axis(1) pstyle(p1)) ///
(epa, keep(1:) ci(ci) ciopts(col(maroon%30) lw(0)) axis(1) pstyle(p2)) ///
, at(at) recast(line) cirecast(rarea) graphreg(c(gs15) m(t+3 l-3 r-0)) ///
plotreg(m(l+0 r+1 t+15)) ytit("density", m(l+1 r+0)) name(d`var', replace) ///
ysc(range(0 1) noex) ylab(, labs(*1) angle(none) nogrid) ytit("Density", m(r-1) s(*1.2)) ///
xsc(range(1 5) noex) xmt(1.1(.1)4.9) xlabel(1(.5)5, nogrid) /// legend(off) ///
xtit("`vart'", m(t+0 b-2) s(*1.2)) legend(order(4 "females" 2 "males") ring(0) pos(11) region(lc(none)) symx(*.8)) ///
xline(`Llsd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`L2sd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`L3sd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`L4sd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`Ulsd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`U2sd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`U3sd' , lc(gs13) lw(*1) lp(shortdash)) ///
xline(`pL' , lc(gs10) lw(*1) lp(shortdash)) /// S-index lower bound (selected)
xline(`pM' , lc(gs10) lw(*1) lp(dash)) /// S-index population mean
xline(`pU' , lc(gs10) lw(*1) lp(shortdash)) /// S-index upper bound (selected)
text(1.07 `pM' "M") ///
text(1.07 `L4sd' "-4 SD") text(1.07 `L3sd' "-3 SD") text(1.07 `L2sd' "-2 SD") text(1.07 `Llsd' "-1 SD") ///
text(1.07 `U3sd' "3 SD") text(1.07 `U2sd' "2 SD") text(1.07 `Ulsd' "1 SD")
addplot: (line patFe patx_e if patx_e>=`L3sd', yaxis(2) xaxis(2) lc(stblue)) (line patMe patx_e if patx_e>=`L3sd', yaxis(2) xaxis(2) lc(stred)), ysc(range(-100 70) alt axis(2) noex) xsc(range(1 5) alt axis(2) noex off) ytit(`change', orient(horizontal) m(l-6 r-2) s(*1.2) axis(2)) xtit("", axis(2)) legend(order(4 "females" 2 "males") ring(0) pos(11) region(lc(none)) symx(*.8))
which results in:
The difference between the two plots is in this case not that large but I can think of several examples using my own data that does.
Ben, it is nice to be able to replicate the figure(s) that support the proposed S-index but I am at a loss what to do to compute it using the published formula:
Maybe I am wrong, but, time permitting, I think that it would not take too much effort to use his R-code and replicate it as a Stata command which would possibly offer more flexibility and functionality, or embed it as an option in one of your own packages (like dstat or cdist).
For example, certainly, I would like to be able to use survey weigths while calculating this S-index.
One methodological question I have about the method and formulation.
Marco strives to extend the ability to compare groups by their averages with the S-index by including (using) the proportion of group members at the extremes of the distribution as set by a certain standard deviation (some substantive choice the analyst should make).
The appeal of the S-index is that it ranges between -1 and 1.
But, if I understand him correctly, would it be possible to extend the above formula to include multiple points along the distribution(?).
Like calculating at each decile or at the center point of the 9 stanines?
Possibly we would then get some sort of Distribution-index that maybe would allow to compare the distribution between two groups also over not-normally distributed data.
Well, if you think that I am 'lost in the woods' then also I am glad to hear from you what you think of the above.
Best regards,
Eric
Dear Ben,
I did some further work on the replication of one of the figures published by:
Del Giudice, M. (2023). The S-index: Summarizing patterns of sex differences at the distribution extremes. Personality and Individual Differences, 205, 112088. https://doi.org/10.1016/j.paid.2023.112088
The objective of the so-called S-index, proposed by Marco Del Giudice, is 'to compare groups that are about equally represented in the population, which makes it ideally suited to the analysis of [on] sex differences' but not limited to that.
You are in a better position than me to assess the mathematical and statistical appropriateness of what Marco proposes.
But, for my research, which does meet the relevant assumptions that he posits, it might be a useful addition to the set of tools for distributional analysis and comparison between groups, like females and males.
So, using
dstat, I have worked on his example data and can now replicate one of the traits of Figure 1 (Openness).I suppose that my code is something of a 'hack' that you might expect from a not too experienced coder like me.
Therefore, to replicate his figure:
using
coefplot(somewhat improving it in the process) I go through the following:The paper and the dta file that I created using his R-code and data, I provide here in a zipped file:
PID 2023 Guidice - R dta.zip
Next, I use this code:
which results in:
My opinion is that this plot is more clear, because of the visualization of the two distribution curves separate from the two 'patterns of proportions' with their scales. Also more informative is that
dstatallows for the visualization of confidence intervals which allows for a further analysis of the distributional differences between the two groups.Also note that Marco applies the Gaussian distribution function which, I suppose, follows from his assumption that the scale data under inspection more or less originates from an 'idealized' normal distribution. That might be the case, or not, and I also (tried to) advance his result but using the Epanechnikov distribution function from
dstat, like:which results in:
The difference between the two plots is in this case not that large but I can think of several examples using my own data that does.
Ben, it is nice to be able to replicate the figure(s) that support the proposed S-index but I am at a loss what to do to compute it using the published formula:
Maybe I am wrong, but, time permitting, I think that it would not take too much effort to use his R-code and replicate it as a Stata command which would possibly offer more flexibility and functionality, or embed it as an option in one of your own packages (like
dstatorcdist).For example, certainly, I would like to be able to use survey weigths while calculating this S-index.
One methodological question I have about the method and formulation.
Marco strives to extend the ability to compare groups by their averages with the S-index by including (using) the proportion of group members at the extremes of the distribution as set by a certain standard deviation (some substantive choice the analyst should make).
The appeal of the S-index is that it ranges between -1 and 1.
But, if I understand him correctly, would it be possible to extend the above formula to include multiple points along the distribution(?).
Like calculating at each decile or at the center point of the 9 stanines?
Possibly we would then get some sort of Distribution-index that maybe would allow to compare the distribution between two groups also over not-normally distributed data.
Well, if you think that I am 'lost in the woods' then also I am glad to hear from you what you think of the above.
Best regards,
Eric