Geometric means - gmci

[Sabrina Gravel]

Hi! I use the command gmci often because I work with a lot of log-normally distributed datasets.

The option by() is not allowed, and I am wondering if it would be easy to add it in the program. I am generally not comfortable with programming but I know that sometimes it's not that complicated to add functions.

I know that ameans allow it, but I quite like gmci!

[Nick Cox]

gmci is a community-contributed command, as you are asked to explain (FAQ Advice #12). I didn't recognise it but search finds a 1998 publication in https://www.stata.com/products/stb/journals/stb41.pdf in which the authors acknowledge earlier work by William Gould and myself.

I don't know how likely the authors are to want to revisit their program 23 years later, but it seems that a minimal hack is to change the code for saved results and then use the machinery of statsby.

Here's a sample:

. sysuse auto, clear
(1978 automobile data)

* original program yields these results
. gmci mpg
 
  Variable  |    Obs     Geometric Mean           [95% Conf. Interval]
------------+-----------------------------------------------------------
       mpg  |     74          20.58444          19.38034      21.86335
 

* gmci2 is a hack

. statsby , by(foreign) : gmci2 mpg
(running gmci2 on estimation sample)

      Command: gmci2 mpg
           ub: r(ub)
           lb: r(lb)
        gmean: r(gmean)
            N: r(N)
           By: foreign

Statsby groups
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..

. l

     +------------------------------------------------+
     |  foreign         ub         lb      gmean    N |
     |------------------------------------------------|
  1. | Domestic   20.59307   18.09168   19.30189   52 |
  2. |  Foreign    26.9354   21.32215   23.96499   22 |
     +------------------------------------------------+

Here's

gmci2

*! 2.0.0 NJC hack 29 June 2021
*! version 1.0.3 SV/JBC 17 Sept 1997 after Ramalheira/Gould/Cox STB-41 sg75
program define gmci2, rclass
    version 5.0
    local varlist "req ex"
    local if "opt"
    local in "opt"
        local weight "aweight fweight"
    local options "Level(int $S_level) Add(real 0) Only"
    parse "`*'"

    parse "`varlist'", parse(" ")
    if `add' < 0 {
    di in r "A negative value cannot be added."
    exit 198
    }
    if `level'<10 | `level'>99 {
      dis in red "level() invalid"
      exit 198
      }
    tempvar var1 junk

    local i 0
    while "`*'"!="" {
        qui gen `var1' = `1' + `add'
        qui su `var1' [`weight' `exp'] `if' `in'
        local min=_result(5)
          if `min' <= 0 {
           if "`if'"=="" {
            qui su `1' [`weight' `exp'] `in' if `1'>0
            }
           else {
            qui su `1' [`weight' `exp'] `if' & `1'>0 `in'
            }
        local minpos=_result(5)
          if `add' > 0 {
                di in r "Nonpositive values still encountered in variable `1'."
            di in r "Minimum value of (`1'+`add') is " %4.1g `min' /*
             */", minimum positive value of `1' is " %4.1g `minpos'
                  }
          if `add' == 0 {
                di in r "Nonpositive values encountered in variable `1'."
            di in r "Minimum value of `1' is " %4.1g `min' /*
             */", minimum positive value of `1' is " %4.1g `minpos'
                  }
        local i 1
            }
        drop `var1'
        mac shift
    }

    if `i' == 1 {
        exit 411
      }

    parse "`varlist'", parse(" ")

    di " "
    di in gr "  Variable  |    Obs     Geometric Mean           [" /*
        */ `level' "%"   /*
        */ " Conf. Interval]"
    di in gr "------------+" _dup(59) "-"

    local msg 0

    while "`1'" ~= "" {
    qui su `1' [`weight' `exp'] `if' `in'
      if _result(5) > 0 & "`only'"~="" {
        qui gen `var1' = `1'
        local star  ""
                }
      else {    
           qui gen `var1' = `1' + `add'
        local msg 1
        local star *
            }
        qui gen `junk' = ln(`var1') `if' `in'
        qui ci `junk' [`weight' `exp'] `if' `in' , /*
            */ `options' level(`level')
        if `add'==0 { local star = "" }
        local namelen = 10 - length("`1'")
        display in ye _dup(`namelen') " " "`1'" _col(13) /*
        */in gr "|"  in ye _col(16) %5.0f _result(1)  /*
        */ in ye _col(27) %12.5f exp($S_3) /*
        */ in ye "      " in ye _col(45) %12.5f exp($S_5) /*
        */ in ye _col(59) %12.5f exp($S_6) in w " `star'"

    return scalar N = _result(1)
    return scalar gmean = exp($S_3)
        global S_4 = exp($S_4)
        return scalar lb = exp($S_5)
        return scalar ub = exp($S_6)
    drop `var1' `junk'
    mac shift
    }

di " "
if `msg'~=0 & `add' >0 {
    di in gr "(" in w "*" in gr ") `add' was added to the " /*
        */ "variable(s) prior to calculating the results"
    }
    
end

[Sabrina Gravel]

Thank you so much!
This allowed me to explore Statsby, which I don't use often. However using this approach I was unable to use more than one variable at once.
I found an easy workaround, which I admittedly should have thought of before posting on the forum.

. sysuse auto, clear
(1978 Automobile Data)

. forval i = 0/1 {
      gmci mpg weight if foreign==`i'
           }
 
  Variable  |    Obs     Geometric Mean           [95% Conf. Interval]
------------+-----------------------------------------------------------
       mpg  |     52          19.30189          18.09168      20.59306 
    weight  |     52        3239.02313        3040.53616    3450.46737 
 
 
  Variable  |    Obs     Geometric Mean           [95% Conf. Interval]
------------+-----------------------------------------------------------
       mpg  |     22          23.96499          21.32215      26.93540 
    weight  |     22        2280.79831        2110.24344    2465.13783

I only wish that I was able to display the values (better yet - the labels) of `i' before each table, as some of my categorical variables have many categories.

[Nick Cox]

It could be a project for you to hack further at the code: there is no incentive like really wanting a command yourself.

However, as the official command ameans is documented to supersede gmci, that would be your own choice.

Those glancing at the manual entry for ameans will see a reference to a 1911 paper by J.M. Keynes, who is better known for some other work in economics.

[Attaullah Shah]

You can also try asrol (from SSC).

ssc install asrol

sysuse auto, clear
bys foreign : asrol mpg, stat(gmean)

list mpg_gmean in 1

     +-----------+
     | mpg_gmean |
     |-----------|
  1. | 19.301893 |
     +-----------+

. list mpg_gmean in 74

     +-----------+
     | mpg_gmean |
     |-----------|
 74. | 23.964987 |
     +-----------+

* Confirm using ameans command

ameans mpg if foreign == 0
    Variable |    Type             Obs        Mean       [95% Conf. Interval]
-------------+---------------------------------------------------------------
         mpg | Arithmetic           52    19.82692        18.50638   21.14747
             |  Geometric           52    19.30189        18.09168   20.59306
             |   Harmonic           52    18.80351        17.67092   20.09123
-----------------------------------------------------------------------------

. ameans mpg if foreign == 1

    Variable |    Type             Obs        Mean       [95% Conf. Interval]
-------------+---------------------------------------------------------------
         mpg | Arithmetic           22    24.77273        21.84149   27.70396
             |  Geometric           22    23.96499        21.32215    26.9354
             |   Harmonic           22    23.18653        20.74585   26.27807
-----------------------------------------------------------------------------

asrol can find the following statistics in rolling windows and by groups (your example involves groups )

 
 -------------------------------------------------------------------------------- 
       sd         Estimates the standard deviation of non-missing values

       mean       Finds the arithmetic mean of non-missing values

       gmean      Finds the geometric mean of positive values

       sum        Adds all the numbers in a given window

       product    Multiplies all the numbers in a given window

       median     Returns median of non-missing values

       count      Counts the number of non-missing observations in a given window

       missing    Counts the number of missing values in a given window

       min        Returns the smallest value in a given window

       max        Returns the largest value in a given window

       first      Returns the first observation in a given window

       last       Returns the last observation in a given window

       perc(k)    Returns the k-th percentile of values in a range. This option
                    must be used in combination with the option median.
--------------------------------------------------------------------------------

asrol is blinking fast compared to rolling command. More details related to asrol can be found here https://fintechprofessor.com/2017/10...tics-in-stata/

[Nick Cox]

Here's a way to show headings.

sysuse auto, clear

forval i = 0/1 {
      local which: label (foreign) `i'
      di "{title:`which'}"
      gmci mpg weight if foreign==`i'
      di 
}

[Sabrina Gravel]

Nick - True! Creating (or hacking) programs for myself is definitely on my "to learn" list. Thank you for the display command in the loop, this is perfect.

Attaullah - Thank you, I will explore this command further.


Working my way up this proverbial learning curve...