Age at Menarche: Warsaw, Poland 1965

Documentation here

Human Biology, Vol. 38, No. 3 (September, 1966), pp. 199-203

Data

The data from this article are included in the MASS package in `R’

ds=menarche; head(ds); tail(ds)

##     Age Total Menarche
## 1  9.21   376        0
## 2 10.21   200        0
## 3 10.58    93        0
## 4 10.83   120        2
## 5 11.08    90        2
## 6 11.33    88        5

##      Age Total Menarche
## 20 14.83   102       95
## 21 15.08   122      117
## 22 15.33   111      107
## 23 15.58    94       92
## 24 15.83   114      112
## 25 17.58  1049     1049

Fit a Normal (GAUSSIAN) Distribution, via probit link in a GLM

probit.fit = 
   glm(cbind(Menarche, Total - Menarche) ~ Age,
         family=binomial(link = probit), 
         data = menarche )

print(probit.fit)

## 
## Call:  glm(formula = cbind(Menarche, Total - Menarche) ~ Age, family = binomial(link = probit), 
##     data = menarche)
## 
## Coefficients:
## (Intercept)          Age  
##    -11.8189       0.9078  
## 
## Degrees of Freedom: 24 Total (i.e. Null);  23 Residual
## Null Deviance:       3694 
## Residual Deviance: 22.89     AIC: 110.9

—- Find 50:50 age (median and mean), i.e. the age where probit = Z = 0, as well as SD

\(\pi[a]\) = Prob[already reached menarche by age a] = \(\Phi[ \frac{age \ - \ \mu}{\sigma} ] = \Phi[ \frac{1}{\sigma} age \ - \ \frac{\mu}{\sigma}]\), so \[\Phi^{-1}[ \ \pi[a] \ ] = \frac{1}{\sigma} age - \frac{\mu}{\sigma} = \beta_1 \ age + \beta_0,\] where \(\beta_0 = -\mu / \sigma\) and \(\beta_1 = 1/\sigma,\) so that \(\sigma = 1/\beta_1\) and \(\mu = - \beta_0/\beta_1\).

Would report mean and sd. to at most 1 decimal place, but using 2 decimal places here for comparison with logit fit

age.50 = -probit.fit$coefficients[1]/
          probit.fit$coefficients[2]

sd = 1/probit.fit$coefficients[2] 


print(noquote(paste("Fitted Mean/Median (SD): ",
            sprintf("%3.2f",age.50),
            " (",
            sprintf("%3.2f",sd),
            ") years",sep="") ))

## [1] Fitted Mean/Median (SD): 13.02 (1.10) years

Fit a Logistic Curve / Distribution, via logit link in a GLM

logit.fit = 
   glm(cbind(Menarche, Total - Menarche) ~ Age,
         family=binomial(link = logit), 
         data = menarche )
logit.fit

## 
## Call:  glm(formula = cbind(Menarche, Total - Menarche) ~ Age, family = binomial(link = logit), 
##     data = menarche)
## 
## Coefficients:
## (Intercept)          Age  
##     -21.226        1.632  
## 
## Degrees of Freedom: 24 Total (i.e. Null);  23 Residual
## Null Deviance:       3694 
## Residual Deviance: 26.7  AIC: 114.8

—- Find 50:50 age (median and mean), i.e. the age where logit = 0, as well as SD

\[\pi[a] = \frac{e ^{\beta_0 + \beta_1 age} }{1+ e ^{\beta_0 + \beta_1 age}} = CDF_{logisticDistrn.}[age]\]

**See Logistic DISTRIBUTION (Wiki)

age.50 = -logit.fit$coefficients[1]/
          logit.fit$coefficients[2]   ; 

sd = ( pi/sqrt(3) ) /logit.fit$coefficients[2] ; 

print(noquote(paste("Fitted Mean/Median (SD): ",
            sprintf("%3.2f",age.50),
            " (",
            sprintf("%3.2f",sd),
            ") years",sep="") ))

## [1] Fitted Mean/Median (SD): 13.01 (1.11) years

Age.at.Menarche = seq(8,18,.1)

y=dlogis(Age.at.Menarche,age.50,1/logit.fit$coefficients[2])
plot(Age.at.Menarche,y,type="l",ylab="Density",
     main="[G]AUSSIAN and [L]OGISTIC DISTRIBUTIONS",
     cex.lab=1.5)
points(Age.at.Menarche,y,pch="L",cex=0.8)

y = dnorm(Age.at.Menarche,age.50,sd)
lines( Age.at.Menarche,y,col="red" )
points(Age.at.Menarche,y,col="red",pch="G",cex=0.8 )

y=plogis(Age.at.Menarche,age.50,1/logit.fit$coefficients[2])
plot(Age.at.Menarche,y,type="l",ylab="CDF",cex.lab=1.5,
     main="[G]AUSSIAN and [L]OGISTIC CUMULATIVE DISTRN. FUNCTIONS")
points(Age.at.Menarche,y,pch="L",cex=0.8)

y = pnorm(Age.at.Menarche,age.50,sd)
lines( Age.at.Menarche,y,col="red" )
points(Age.at.Menarche,y,col="red",pch="G",cex=0.8 )
abline(h=seq(0.0,1.0,0.1)); abline(h=0.5,lwd=2)
abline(v=seq(9,17,1))