fork download
copy 
  1. # Must have R, RStudio, RTools, and JAGS installed.  I have R3.3.1, RSTudio 0.97.551, RTools 3.4, and JAGS 4.2.0.
  2. # Must install.package(rjags)
  3. # 
  4. # Model
  5. # K = number of species in sample
  6. # N = number of possible species in the sample
  7. # M = number of short reads
  8. # theta[N] = Indicator of species in sampe, species are numbered 1-N; 0 out, >0 in; prior multinomial(priorspecies,K)
  9. # alpha[N]= Abundance of species in the sample; 
  10. # prior dirichlet with those in sample have prior proportion 1/K thse not in sample have prior proportion 1/N
  11. # z[M] = latent, unknown assignment of reads to species; possible values are 1-N; each read is assigned a value using alpha[] probabilities
  12. # r[M] = list of BLAST alignments associated with each reads
  13. # r[m] ~ multinomial(length(r[m], coMatch[z[m],]); this is currently one read and a trick is used to get the likelihood
  14. # coMatch[N, N] = given true species (row), the likelihood of aligning to another species (col)
  15.  
  16.  
  17. library(rjags)
  18.  
  19. K <- 5
  20. N <- 100
  21. M <- 1000
  22. theta <- rep(0,N)
  23. theta[1:K] <- 1
  24. alpha <- rep(0,N)
  25. alpha[theta] <- 1/K
  26.  
  27. # First example has only one BLAST return that includes 
  28. # 1 (1/50th), 
  29. # 2 (9/50th), 
  30. # 3 (2/50th), 
  31. # 4 (2/50th), 
  32. # 5 (2/50th), 
  33. # 6 (10/50th), 
  34. # 8 (2/50th), 
  35. # 9 (3/50th), 
  36. # 10 (8/50th), and 
  37. # 50 (11/50th)
  38. r <- sample(c(1,2,2,2,2,2,2,2,2,2,3,3,4,4,5,5,6,6,6,6,6,6,6,6,6,6,8,8,9,9,9,10,10,10,10,10,10,10,10,50,50,50,50,50,50,50,50,50,50,50), size=M, replace=T)
  39.  
  40. # Sepcify a simple coMatch matrix that makes all off-diagonal entries low probability
  41. likelihood <- matrix(0.1,nrow=100, ncol=100)
  42. diag(likelihood) <- rep(1,100)
  43.  
  44.  
  45. write('
  46.       model {
  47.       # Likelihood:
  48.       for( i in 1 : M ) {
  49.         zeros[i] ~ dpois(phi[i])         # likelihood is exp( - phi[i])
  50.         phi[i] <- -1*log(likelihood[r[i], z[i]])         
  51.         z[i] ~ dcat(alpha[])  
  52.       }
  53.  
  54.     # Prior:
  55.  
  56.     alpha ~ ddirch(abundance)
  57.     for(j in 1:N){
  58.       abundance[j] <- min(theta[j]*1/K, 1/N)
  59.       }
  60.     theta ~ dmulti(priorSpecies, K)
  61.  
  62.     }', "mixturemodel4.bug")
  63.  
  64. # Hyperparameters  	  
  65. #  priorAbundance <- rep(1/N,N)
  66. #  priorAbundance[theta] <- 1
  67.   priorSpecies <- rep(1/N,N)
  68.  
  69. # false data variable for implementing a vector likelihood function in rjags
  70.   zeros <- rep(0,M)
  71.  
  72.  
  73. # Create a data list
  74. dataList <- list(
  75.   "r"= r,
  76.   "zeros" = zeros,
  77.   "likelihood" = likelihood,
  78.   "K" = K,
  79.   "N" = N,
  80.   "M"= M, 
  81.   "priorSpecies" = priorSpecies)
  82.  
  83. # List of parameters to be monitored  
  84. parameters <- c(
  85.   "alpha",
  86.   "theta")
  87.  
  88. # Set initial values
  89. initsValues <- list(
  90.   "alpha" = alpha,
  91.   "theta" = theta)
  92.  
  93.  
  94. # JAGS Set-up
  95. adaptSteps <- 50              #number of steps to "tune" the samplers
  96. burnInSteps <- 50             #number of steps to "burn-in" the samplers
  97. nChains <- 1                    #number of chains to run
  98. numSavedSteps <- 50           #total number of steps in chains to save
  99. thinSteps <- 1                  #number of steps to "thin" (1=keep every step)
  100. nIter <- ceiling((numSavedSteps*thinSteps )/nChains) 	#steps per chain
  101.  
  102. # Create, initialize, and adapt the model
  103. jagsModel <- jags.model("mixturemodel4.bug", 
  104.                         data=dataList, 
  105.                         inits=initsValues, 
  106.                         n.chains=nChains, 
  107.                         n.adapt=adaptSteps)
Success #stdin #stdout #stderr 0.22s 60752KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
stderr
copy 
Error in library(rjags) : there is no package called ‘rjags’
Execution halted
https://ideone.com/X6TlGM
language:
R (R 3.5.2)
created:
7 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!