rm(list=ls(all=T)) set.seed(3) library("topicmodels") setwd('I:/Forestry Users/drvalle/lda/ecol letters/rcode') dat.agg=data.matrix(read.csv('fake data agg.csv',as.is=T)) #sampling unit x species matrix #run LDA SEED=2010 VEM=LDA(dat.agg,k=3, control = list(seed = SEED)) #get parameter estimates z=posterior(VEM) commun.plot=z$topics commun.spp=z$terms #plot relative abundance of component communities for each sampling unit 1,...,1000 (i.e., along the gradient) par(mfrow=c(1,1)) plot(NA,NA,xlim=c(0,1000),ylim=c(0,1),xlab='Gradient/Sampling units',ylab='Relative abundance') for (i in 1:3){ lines(1:1000,commun.plot[,i],col=i) } #plot relative abundance of species 1,...,2000 in component community par(mfrow=c(3,1)) for (i in 1:3){ plot(1:200,commun.spp[i,],type='l',col=i,ylim=c(0,0.02),xlab='Species',ylab='Relative abundance') }