Structural model
- Maturity is a funciton of age, \(p(a)\)
- \(p(a)\) should be a monotonically increasing function of \(a\) like a sigmoid 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 \]
\[ p(a)=\displaystyle \int_{-\infty}^a f(z)dz \]
# 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")
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")
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")
# 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)