Maturity model

Structural model
- Maturity is a funciton of age, \(p(a)\)
- \(p(a)\) should be a monotonically increasing function of \(a\) like a sigmoid function:

  1. Logistic function

\[ p(a)=\displaystyle \frac{e^{\beta_0 + \beta_1 a}}{1+e^{\beta_0 + \beta_1 a}} \ \Longrightarrow \ \log \frac{p(a)}{1-p(a)} = \beta_0 + \beta_1 a \]

  1. Probit function

\[ p(a)=\displaystyle \int_{-\infty}^a f(z)dz \]

Drawing logistic function

# Definition of logistic function (the parameters are beta0 and beta1)
LogisticFn <- function(x,beta0,beta1) {1.0/(1.0+exp(-beta0-beta1*x))}
Age <- seq(0,10,0.01)

# beta0=-5, beta1=1
plot(Age, LogisticFn(Age,-5,1), type="l", ylab="Maturity", lty=1)

# beta0=-2.5, beta1=0.5
points(Age, LogisticFn(Age,-2.5,0.5), type="l", lty=2)

# beta0=-10, beta1=2
points(Age, LogisticFn(Age,-10,2), type="l", lty=3)

# legends
legend(0, 1, lty=1:3, 
 legend=c("beta0=-5, beta1=1","beta0=-2.5, beta1=0.5","beta0=-10, beta1=2"),bty="n")

Statistical model and generation of example data

Distributional assumption
- Two state (mature/immmature) -> Binomial distribution
- \(N_a\): number of sampled individuals of age \(a\)
- \(Y_a\): number of mature individuals among \(N_a\) individuals
- \(p(a)\): maturity (probability) in age \(a\)
- \(Y_a \sim Bin(N_a, p(a))\)

set.seed <- 2018

Amin <- 1
Amax <- 10
AgeRange <- Age <- Amin:Amax
NAgeClass<- length(AgeRange)

# True beta0=-5, beta1=1
Nvec <- rep(10, NAgeClass)
pvec <- LogisticFn(AgeRange, -5, 1)
pvec
##  [1] 0.01798621 0.04742587 0.11920292 0.26894142 0.50000000 0.73105858
##  [7] 0.88079708 0.95257413 0.98201379 0.99330715
# Generation of simulation data
Yvec <- rbinom(NAgeClass, Nvec, pvec)
Yvec
##  [1]  0  1  2  1  4  6  7  9 10 10
Maturity.obs <- Yvec/Nvec
Maturity.obs
##  [1] 0.0 0.1 0.2 0.1 0.4 0.6 0.7 0.9 1.0 1.0
# Plotting the data
plot(Age, LogisticFn(Age,-5,1), type="l", ylab="Maturity")
points(Age, Maturity.obs, pch="*")
legend(6, 0.4, lty=1:0, pch=c("","*"), legend=c("True","Observed"), bty="n")

Statistical estimation of maturity curve using “glm” function

Res.glm.logit <- glm(cbind(Yvec,Nvec-Yvec)~Age, family=binomial)
Res.glm.logit
## 
## Call:  glm(formula = cbind(Yvec, Nvec - Yvec) ~ Age, family = binomial)
## 
## Coefficients:
## (Intercept)          Age  
##      -4.697        0.854  
## 
## Degrees of Freedom: 9 Total (i.e. Null);  8 Residual
## Null Deviance:       69.98 
## Residual Deviance: 4.395     AIC: 24.65
# Summary
summary(Res.glm.logit)
## 
## Call:
## glm(formula = cbind(Yvec, Nvec - Yvec) ~ Age, family = binomial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.98108  -0.46507   0.04643   0.67068   0.99109  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -4.6970     0.9169  -5.122 3.02e-07 ***
## Age           0.8540     0.1578   5.413 6.21e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 69.9787  on 9  degrees of freedom
## Residual deviance:  4.3952  on 8  degrees of freedom
## AIC: 24.654
## 
## Number of Fisher Scoring iterations: 4
# Parameter estimates
Res.glm.logit$coef
## (Intercept)         Age 
##  -4.6969770   0.8539958
beta0.est <- Res.glm.logit$coef[1]
beta1.est <- Res.glm.logit$coef[2]

# CI of the paramters
confint(Res.glm.logit)
## Waiting for profiling to be done...
##                  2.5 %    97.5 %
## (Intercept) -6.7455369 -3.101847
## Age          0.5805778  1.207944
# Visual presentation of results
plot(Age, LogisticFn(Age,-5,1), type="l", ylab="Maturity")
points(Age, Maturity.obs)
points(Age, LogisticFn(Age,beta0.est,beta1.est), type="l", col="red")

Confidence interval for the curve

# Estimates
pred.se<-predict (Res.glm.logit, se.fit=TRUE)  #Linear predictor

#Linear predictor's CI
Lower <- pred.se$fit - qnorm(0.975, 0, 1)* pred.se$se.fit
Upper <- pred.se$fit + qnorm(0.975, 0, 1)* pred.se$se.fit

#Logistic transformation of CI
logit.t<-function(x){ 1/(1+exp(-x)) }
pred.logit<- sapply(pred.se$fit, logit.t)
Lower.logit <- sapply(Lower, logit.t)
Upper.logit<- sapply(Upper, logit.t)
pcon.logit<-cbind(pred.logit, Lower.logit, Upper.logit)
plot(AgeRange, Maturity.obs, xlim=c(0,10), ylim=c(0,1))
matlines(AgeRange, pcon.logit, lty=c(1,2,2), col="blue", lwd=1.5)