Which GLM to use for analyzing the trend in the sex ratio between different years? I have data on sex ratios (number of males and the number of females) for 13 years. The samples are independent because they are hunted individuals.
I want to see if there is a trend of decreasing or increasing males or females over the years.
I did a Chi-square test and now I know that the sex ratios differ between years. By plotting the sex ratios I know there is a trend that could be explained by a quadratic term (females increase at the beginning and then they decrease), however, I would like to use a GLM to test this.
What is the best GLM for this?
year <- c(2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 
           2012, 2014, 2015, 2016)

females <- c(4503, 4500, 3345, 3408, 1903, 4398, 2784, 4607, 3714, 
              3934, 1488, 551, 1591, 2391)

males <- c(1359, 1360, 916, 827, 561, 745, 537, 759, 367, 492, 682, 
            112, 455, 959)

 A: The Beta is not the appropriate way to model this in my opinion.  Using a GAM.  Let's set up the problem
    library(tidyverse)
    library(mgcv)
    year <- c(2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 
              2011,  2012, 2014, 2015, 2016)
    
    females <- c(4503, 4500, 3345, 3408, 1903, 4398, 2784, 4607, 
                 3714, 3934, 1488, 551, 1591, 2391)
    
    males <- c(1359, 1360, 916, 827, 561, 745, 537, 759, 367, 492, 
               682, 112, 455, 959)
    
    N = females + males
    
    d = tibble(year = year - min(year), y = females, N = N)
        
    model = gam(cbind(y, N-y) ~ s(year), data = d, 
                family = binomial())

Then, we can plot the model on the probability scale by calling plot.  If you just call plot(model), the model will plot on the scale of the linear predictor.  Since we are using a binomial model, you can call plot(model, trans = plogis) to plot on the scale of the probability.  If you call this, you get

The plot here shows the proportion of females in the population (since I set y to the number of females). This trend looks fairly non-linear to me, and I don't think a quadratic term properly captures this.  I went ahead and fit a glm anyway by doing
    model2 = glm(cbind(y, N-y) ~ poly(year,2), data = d, 
                 family = binomial())

and then computed the AIC for each model.  The model for the GAM is much smaller, (240.01 as compared to 737.37) so between the two models the GAM seems to take the cake.
Is this helpful?
