I am conducting a topic modeling study via Latent Dirichlet Allocation (LDA) on scientifc abstracts over a range of about 20 years (in R). One of the outputs of LDA is a document-topic-distribution matrix, columns represent topics, rows show the probability of each topic per document. Each entry in a row may be a value between 0 and 1, while the probabilities sum up to 1 in each row. Using the information on the publication date of the documents might be used to derive timely trends of individual topics. Others, e.g., Griffiths and Steyvers (2004), have simply used the mean of the topic probabilities per year to derive a trend. It should be clear that the expected probability values per topic per year are very small values, since several of the entries for an individual topic and a given year are close to zero.

My question now is: How might a Generalized Additive Model (GAM) best be used to create a smoothed trend (rather for "simple" visualization of trends than for prediction)?

Using a GAM instead of mean values makes sense from my understanding (please correct me if I am wrong), since usually the number of documents available increases with time, so that there are many more data points in recent than in past years - using mean values to build trends somehow discards this information and the steps such as minimizing erros, etc.

For an LDA model with 300 topics for 26533 documents I have fitted GAMs for each topic and I am not sure how to deal with the quality of the results. Maybe someone can support me in improving my understanding. I am relatively new to GAMs so please excuse if I have forgotten to report anything, etc. or maybe have even made wrong choices for the model (just started to read and trying to understand S. Wood's book on GAMs with R).

In the following a model summary for one selected topic. The resulting explained deviance is very low, however, the GCV is very low and the coefficients seem to be significant (see model summary below). This seems a bit contradictory, but probably I have not fully understood the values, yet.

I would appreciate any hints, explanations, criticism.

UPDATE: In response to the answer of @Gavin Simpson I have fitted a gam with family = betar and provided some more details on the data. I have further included a link to the data of the topic used for the models in below code.

I think I have not fully understood how the different models are related regarding prediction of the response. The betar based models appear as a straight line at the scale of the predictions of the other gaussian based models. However, when plotting the former models at their "own" scale, the course of the curve makes sense, especially considering the number of documents per year, making data before 2005 less relevant during fitting (at least, this is my understanding). I am more lost than before knowing about the use of the beta distribution...

Regarding Dirichlet Regression I am still trying to work my way through and might provide an additional update.

# the corresponding file can be found here
# https://github.com/manuelbickel/textility/raw/master/tests/topic_data.RDS
topic_data = readRDS("topic_data.RDS")
# for a first impression plot a histrogram /  the majority of values are zero...
# hence we are some kind of interested in the outliers
# some numbers
#      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# 0.000000 0.000000 0.000000 0.004318 0.000000 0.487500 

enter image description here

# plot models, etc. - the final combined plot is at the end
# plot the mean values
# NOTE: the models will not be fit on these mean values but using all data points
plot(topic_data[, mean(probability), by = year])
# rough indication to know which years contain the most information,i.e., number of documents per year
# of course, we have to scale the number down to the scale of the probability
# the number of docs increase exponentially, that is another challenge for the models, I guess
lines(topic_data[, .N, by = year][, N:=N/7e05], lty = 4)
year_range = data.frame(year = topic_data[, unique(year)])
# check the linear trend
lines(year_range$year, predict(topic_data[, lm(probability ~ year)], newdata = year_range))
# try different gam models
# k will be estimated automatically, just put -1 explicitly to be clear
m1 = topic_data[, mgcv::gam(probability ~ s(year, bs = "cs", k = -1), method="GCV.Cp")]
m2 = topic_data[,mgcv::gam(probability ~ s(year, bs = "cs", k = -1), method="REML")]
m3 = topic_data[,mgcv::gam(probability ~ s(year, bs = "cs", k = -1), family = betar, method="GCV.Cp")]
m4 = topic_data[,mgcv::gam(probability ~ s(year, bs = "cs", k = -1), family = betar, method="REML")]
# add models with gaussian distribution to plot
lines(year_range$year, predict(m1, newdata = year_range), col = 1)
lines(year_range$year, predict(m2, newdata = year_range), col = 2)
# add prediction for on of the models based on beta distribution
lines(year_range$year, predict(m3, newdata = year_range, type = "response"), col = 3)

enter image description here

# plot the models based on beta distribution separately
plot(year_range$year, predict(m3, newdata = year_range, type = "response"), col = 3, type = "p")
lines(year_range$year, predict(m4, newdata = year_range, type = "response"), col = 4) 

enter image description here

# Family: gaussian 
# Link function: identity 
# Formula:
#   probability ~ s(year, bs = "cs", k = -1)
# Parametric coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 0.0043179  0.0001473   29.32   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Approximate significance of smooth terms:
#   edf Ref.df     F p-value    
# s(year) 6.787      9 10.45  <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# R-sq.(adj) =  0.00338   Deviance explained = 0.364%
# GCV = 0.00057552  Scale est. = 0.00057535  n = 26533
# Family: Beta regression(9.622) 
# Link function: logit 
# Formula:
#   probability ~ s(year, bs = "cs", k = -1)
# Parametric coefficients:
#   Estimate Std. Error z value Pr(>|z|)    
# (Intercept) -5.417052   0.006157  -879.8   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Approximate significance of smooth terms:
#   edf Ref.df Chi.sq p-value
# s(year) 0.1044      9   0.01   0.805
# R-sq.(adj) =  -3.52e-06   Deviance explained = 0.0883%
# -REML = -5.8463e+05  Scale est. = 1         n = 26533

1 Answer 1


You are modelling a proportion and you need to fit a model that accounts for the bounded nature of the data.

One option that I have been using recently is to fit a GAM to single topic's proportion using a beta response, where this treats all other topics as complement of the response (i.e. if the response value is 0.3, the remaining 0.7 is from all other topics).

gam(probability ~ s(year, bs = "cs"), data = foo, family = betar)

if you repeat this for all topics, you'll find that the estimated/predicted mean proportions don't add up to 1. For that you'd need to fit a Dirichlet regression, which as far as I can tell can't be fitted via functions in a package available for R whilst retaining the nice features of mgcv to do smoothness selection.

I have used the DirichletReg package to fit the model with bs() or ns() terms from splines package to get smooth trends, but you need to specify the complexity of each smooth.

  • $\begingroup$ Thank you for your valuable explanations. I have updated my answer accordingly and it seems there is a mistake in my code. The family specified in the model does not have an effect. Would you mind to have a look (sorry, I did not intend to get into programming here)? If not, I will still acccept your answer, just wanted to check, if you find time for one more comment/hint/... $\endgroup$ Jul 28, 2018 at 20:52
  • $\begingroup$ And, "Great!", Dirichlet Regression, never heard of it until now, but it sounds perfect for the given context. I have tried to set up the model DirichReg(Y ~ bs(year, df = 9), tt) after using DR_data (df = 9 just for getting started) but it seems I also made a mistake here. I get: "Error in prepareFixed(start = start, activePar = activePar, fixed = fixed) : At least one parameter must not be fixed using argument 'fixed'" Maybe you can also point me to my mistake here, otherwise, I might turn to SO. Many thanks and sorry for bothering you with this. $\endgroup$ Jul 28, 2018 at 20:52
  • $\begingroup$ I realized that my DirichReg model is wrong, but might need some time to understand how it works. Found, e.g., some slides by Marco Maier, have to work though them (and more) first before really being able to use the model. Still, maybe you can point me to some helpful resources / examples. $\endgroup$ Jul 28, 2018 at 21:06
  • $\begingroup$ I found my error for the models with the beta distribution. Simply a typo, sorry. My updated answers leaves me with the question how to interpret the different results. Seems like using betar does not lead to significant smooth terms. Maybe the fact that the earlier years only include very few data points causes this, I am a bit lost, meaning a lot of work/learning ahead... $\endgroup$ Jul 28, 2018 at 22:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.