I have a data set with 40 sites (long*lat) by 2 ages by 2 genders by 6 words with 6 pronunciation categories (coded 0 to 5) as simulated here:

sound         <- sample(0:5, size=960, replace=T)
word          <- as.factor(rep(c('alpha', 'bravo', 'charlie', 'delta', 'echo', 'foxtrot'), times=960/6))
age           <- as.factor(rep(c('young', 'old'), times=960/2))
gender        <- as.factor(rep(c('female', 'female', 'male', 'male'), times=960/4))
long          <- rep(c(runif(40)), each=960/40)
lat           <- rep(c(runif(40)), each=960/40)
pronunciation <- data.frame(sound, word, age, gender, long, lat)

And I want to run a simple multinom GAM (for now), which works fine:

multinom = gam(list(sound ~ word + s(long, lat),
                          ~ word + s(long, lat),
                          ~ word + s(long, lat),
                          ~ word + s(long, lat),
                          ~ word + s(long, lat)),
                          data = pronunciation, family=multinom(K=5))

When I consider the actual distribution of my pronunciation categories:

sound <- sample(0:5, size=960, prob=c(0.17, 0.31, 0.21, 0.28, 0.02, 0.03), replace=T)

Something strange happens. Depending on the simulation run this error pops up:

Warning message:
In newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,  :
  Iteration limit reached without full convergence - check carefully

And with my actual data I get:

Warning message:
In newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,  :
  Fitting terminated with step failure - check results carefully

Is there a possibility to tweak the model in order to run or do I have to modify my data?