Using weights in a multinomial GAM (mgcv)

Question

I'm interested in understanding the factors that affect the proportion of fish that leave a river as fry, parr or smolt (summing to 1). I know its possible to do this using nnet::multinom which works fine, but I wanted to see whether I could fit these relationships with splines using the mgcv package. To do this, I was going to fit a multinomial model with the 3 classes, with regression weights relative to the proportion of that class observed on that day. I can fit a multinomial model with mgcv, however, it doesn't seem to be accepting weights. See below for an example - the outputs are identical with and without weights. Does anyone know whether its possible to use weights in a multinomial gam implemented in mgcv?

library(mgcv)
set.seed(6)
## simulate some data from a three class model
n <- 1000
f1 <- function(x) sin(3*pi*x)*exp(-x)
f2 <- function(x) x^3
f3 <- function(x) .5*exp(-x^2)-.2
f4 <- function(x) 1
x1 <- runif(n);x2 <- runif(n)
eta1 <- 2*(f1(x1) + f2(x2))-.5
eta2 <- 2*(f3(x1) + f4(x2))-1
p <- exp(cbind(0,eta1,eta2))
p <- p/rowSums(p) ## prob. of each category 
cp <- t(apply(p,1,cumsum)) ## cumulative prob.
## simulate multinomial response with these probabilities
## see also ?rmultinom
y <- apply(cp,1,function(x) min(which(x>runif(1))))-1
## plot simulated data...
plot(x1,x2,col=y+3)

## now fit the model...
b <- gam(list(y~s(x1)+s(x2),~s(x1)+s(x2)),family=multinom(K=2))

# Try with weights
reg_weights <- sample(1:100, length(y),replace=T)
b_2 <- gam(list(y~s(x1)+s(x2),~s(x1)+s(x2)),family=multinom(K=2), weights = 
reg_weights)
b; b_2 #identical

Looks like your are trying to model compositional data. Try looking up "[r] estimate dirichlet". I was finding material that suggest that the VGAM package might hold promise. — IRTFM
Brilliant, yes looks like the proportional odds model in vgam does what I want. Thanks! — Hugh Sturrock
You should write up an answer and post it. It's encouraged that questioners post answers to their own questions. — IRTFM

Hugh Sturrock Hugh Sturrock · Accepted Answer · 2018-05-10T23:58:17

Update - I think the VGAM package will do what I want. Here's a worked example using the pneumo data from VGAM to create % in each class.

library(mgcv)
library(VGAM)

# Get pneumo data and log transform exposure time
pneumo <- transform(pneumo, let = log(exposure.time))

# Generate % in each class - note, not suggesting you do 
# this if you have counts, just as an illustrative example if 
# you have % in each class
props <- round(prop.table(as.matrix(pneumo[,c("normal", "mild", 
"severe")]),margin=1)*100,0)

# Fit vgam mod with smoothed relationship with log exposure time
fit <- vgam(props ~ s(let), propodds, data = pneumo)
summary(fit)
plot(fit)

Using weights in a multinomial GAM (mgcv)

1 Answers