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 library(data.table) topic_data = readRDS("topic_data.RDS") # for a first impression plot a histrogram / the majority of values are zero... hist(topic_data$probability) # hence we are some kind of interested in the outliers # some numbers summary(topic_data$probability) # Min. 1st Qu. Median Mean 3rd Qu. Max. # 0.000000 0.000000 0.000000 0.004318 0.000000 0.487500
# 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)
# 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)
summary(m1) # 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 summary(m3) # 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