Finding a reaction norm in R using logistic regression with binomial errors.

Question

I am trying to calculate 'reaction norms' for a fish species. This is essentially the length at which the probability that a fish become mature equals 50% for a particular age class.

I know I have to use a logistic regression model with binomial errors but I can't work out how to calculate this from the summary outputs or plot the regression successfully!

I have a data set that has: 'age' classes in (1,2,3,4,5,6),'Lngth' data in mm and 'Maturity' data (Immature/Mature - 0/1).

I am running a glm as follows

Model<-glm(Maturity~Lgnth, family=binomial(logit)) 

This however does not take into account the different age classes (I would really like to avoid creating whole new data sets for each age classes as I have multiple year ranges to test).

And even so, I do not understand how I interpret the summary output to give me a length at which the probability of being mature equals 50%, along with the standard errors of this figure.

I also can't quite get the code right to plot this. Ideally id have one plot with lngth along the x axis, probability along the y and six lines/curves representing each age classes.

I would really appreciate any help any one could provide! I know this can all be achieved but I am really struggling.

Cheers

Answer 1

I'll demonstrate computing the reaction norms for a two class problem. Let's let $p$ denote the probability of maturity predicted by the model, $L$ the measured length of a fish, and $A$ the class variable (which I'm going to assume is two class, to simplify). Then the logistic model would be something like

$$ \log \left( \frac{p}{1-p} \right) = -1 + .5 L + .5 A $$

where $A$ is a coded indicator, $0$ or $1$, say of "young fish" or "old fish". We want to calculate the cut point at $p = .5$ for both young and old fish. For young, plugging $p=.5$ and $A=0$ into the equation gives:

$$ 0 = -1 + .5 L \Rightarrow L = 2 $$

For old fish $A=1$, so our equation is

$$ 0 = -1 + .5 L + .5 \Rightarrow L = 1 $$

For your plotting question, I'll show you what I learned from Gellman and Hill's book Data Analysis Using Regression and Multilevel/Hierarchical Models. First, let's make some data to model and subsequently plot

N <- 250
expit <- function(t) {(exp(t)/(1+exp(t)))}

X <- data.frame(L=runif(N, 0, 6), A=rbinom(N, 1, .5))
p <- expit(-1 + .5*X$L + .5*X$A)
X$y <- rbinom(N, 1, p)

Fitting a logistic model

model <- glm(y ~ L + factor(A), data=X, family="binomial")
print(model)

Call:  glm(formula = y ~ L + A, family = "binomial", data = X)

Coefficients:
(Intercept)            L     factor(A)1  
    -0.8627       0.4239         0.6917  

If $A$ had more than two classes, you would see additional parameter estimates like factor(A)2 and factor(A)3 etc.

Now here's how to make a nice plot of the data and the probability of maturity for each of the two classes. First, add some jitter to the y values, so they don't crowd and stack over each other

add_jitter <- function(v) ifelse(v == 0, runif(length(v), v, v+.05), runif(length(v), v-.05, v))
y_jitter <- add_jitter(X$y)

Then pull off the coefficients of your model

co <- coef(model)

Then plot

plot(X$L, y_jitter, pch=16, col="grey",
     xlab = "Length of Fish", 
     ylab="Probability of Maturity",
     main="Fishies")
curve(expit(co[1] + co[2]*x), add=TRUE)
curve(expit(co[1] + co[2]*x + co[3]), lty=2, add=TRUE)
legend("right", c("Young Fish", "Old Fish"), lty=c(1, 2), bty="n", lwd=2)

This looks really nice