FMM-2cluster.R

library(readxl)
library(boot)
library(coda)

## Warning: package 'coda' was built under R version 4.0.2

library(R2WinBUGS)

## Warning: package 'R2WinBUGS' was built under R version 4.0.2

# install terlebih dahulu link berikut
# https://sourceforge.net/projects/mcmc-jags/files/
# Untuk Menghitung Nilai DIC

library(R2jags)

## Warning: package 'R2jags' was built under R version 4.0.2

## Loading required package: rjags

## Warning: package 'rjags' was built under R version 4.0.2

## Linked to JAGS 4.3.0

## Loaded modules: basemod,bugs

## 
## Attaching package: 'R2jags'

## The following object is masked from 'package:coda':
## 
##     traceplot

#DATA
dataset1 <- read_excel('Data.xlsx')

var1 <- dataset1[,1:6]
T <- dataset1$TR2

alpha <- c(1,1)
N <- nrow(dataset1)
NLR<-var1$NLR
LMR<-var1$LMR
PLR<-var1$PLR
HPR<-var1$HPR
PWR<-var1$PWR
LWR<-var1$LWR

df<-list("NLR","LMR","PLR","HPR","PWR","LWR","N","alpha","T")

#INIT
inits = function() {
  lambda1 = mean(var1$NLR[1:40]) +rnorm(1,0,.01)
  theta = 60
  sigma2 = var(var1$NLR[1:40])
  return(list(lambda = c(lambda1, NA),
              theta1 = theta,
              tau = 1/sigma2,
              P = c(40, 60-40)/60))
}

#MODEL
mixmodel=function() {
  for( i in 1 : N ) {
    NLR[i] ~ dnorm(mu1[i], tau)
    LMR[i] ~ dnorm(mu2[i], tau)
    PLR[i] ~ dnorm(mu3[i], tau)
    HPR[i] ~ dnorm(mu4[i], tau)
    PWR[i] ~ dnorm(mu5[i], tau)
    LWR[i] ~ dnorm(mu6[i], tau)
    mu1[i] <- lambda[T[i]]
    mu2[i] <- lambda[T[i]]
    mu3[i] <- lambda[T[i]]
    mu4[i] <- lambda[T[i]]
    mu5[i] <- lambda[T[i]]
    mu6[i] <- lambda[T[i]]
    T[i] ~ dcat(P[]) }
  P[1:2] ~ ddirch(alpha[])
  theta1 ~ dnorm(0.0, 1.0E-6)%_%I(0.0, )
  lambda[1] ~ dnorm(0.0, 1.0E-6)
  lambda[2] <- lambda[1] + theta1
  tau ~ dgamma(0.001,0.001)
  sigma <- 1 / sqrt(tau)
}

parameters <- c("lambda","sigma","P")

FMM.fit<-jags(df, inits=inits, parameters, n.chains=1, n.burnin=100000, 
              n.iter=200000, model.file=mixmodel, DIC=TRUE)

## module glm loaded

## Warning in jags.model(model.file, data = data, inits = init.values, n.chains =
## n.chains, : NAs introduced by coercion

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 362
##    Unobserved stochastic nodes: 62
##    Total graph size: 493
## 
## Initializing model

print(FMM.fit)

## Inference for Bugs model at "C:/Users/NALEND~1/AppData/Local/Temp/RtmpYXbxCX/model213c7bb93b48.txt", fit using jags,
##  1 chains, each with 2e+05 iterations (first 1e+05 discarded), n.thin = 100
##  n.sims = 1000 iterations saved
##            mu.vect sd.vect     2.5%      25%      50%      75%    97.5%
## P[1]         0.533   0.242    0.088    0.343    0.537    0.724    0.946
## P[2]         0.467   0.242    0.054    0.276    0.463    0.657    0.912
## lambda[1]   39.977   8.349   19.418   36.147   40.921   45.210   52.428
## lambda[2]   53.563  11.097   38.727   46.934   51.710   57.290   82.213
## sigma       95.919   3.438   89.519   93.448   95.807   98.138  103.272
## deviance  4309.777   2.657 4305.807 4307.924 4309.255 4311.150 4316.209
## 
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 3.5 and DIC = 4313.3
## DIC is an estimate of expected predictive error (lower deviance is better).

# ada notif merah setelah melakukan fitting jangan khawatir langsung PRINT saja
# Koding print(FMM.fit)