lca

lca

library(BayesLCA)

## Warning: package 'BayesLCA' was built under R version 3.6.2

## Loading required package: e1071

## Loading required package: coda

data(Alzheimer)
head(Alzheimer)

##   Hallucination Activity Aggression Agitation Diurnal Affective
## 1             0        0          0         0       0         0
## 2             0        0          0         0       0         0
## 3             0        0          0         0       0         0
## 4             0        0          0         0       0         0
## 5             0        0          0         0       0         0
## 6             0        0          0         0       0         0

fit2 <- blca.em(Alzheimer, 2)

## Restart number 1, logpost = -749.44... 
## Restart number 2, logpost = -749.44... 
## New maximum found... Restart number 3, logpost = -749.42... 
## Restart number 4, logpost = -749.44... 
## Restart number 5, logpost = -749.44...

summary(fit2)

## __________________
## 
## Bayes-LCA
## Diagnostic Summary
## __________________
## 
## Hyper-Parameters: 
## 
##  Item Probabilities:
## 
##  alpha: 
##         Hallucination Activity Aggression Agitation Diurnal Affective
## Group 1             1        1          1         1       1         1
## Group 2             1        1          1         1       1         1
## 
##  beta: 
##         Hallucination Activity Aggression Agitation Diurnal Affective
## Group 1             1        1          1         1       1         1
## Group 2             1        1          1         1       1         1
## 
##  Class Probabilities:
## 
##  delta: 
## Group 1 Group 2 
##       1       1 
## __________________
## 
## Method: EM algorithm  
## 
##  Number of iterations: 24 
## 
##  Log-Posterior Increase at Convergence: 0.0004977783 
## 
##  Log-Posterior: -749.4208 
## 
##  AIC: -1524.842 
## 
##  BIC: -1570.09

fit3<- blca.em(Alzheimer, 3, restarts=20)

## Restart number 1, logpost = -743.98... 
## Restart number 2, logpost = -745.03... 
## Restart number 3, logpost = -745.04... 
## Restart number 4, logpost = -744.28... 
## Restart number 5, logpost = -745.03... 
## Restart number 6, logpost = -744.28... 
## Restart number 7, logpost = -744.68... 
## Restart number 8, logpost = -745.04... 
## Restart number 9, logpost = -745.03... 
## New maximum found... Restart number 10, logpost = -742.79... 
## Restart number 11, logpost = -744.28... 
## Restart number 12, logpost = -744.28... 
## Restart number 13, logpost = -744.29... 
## Restart number 14, logpost = -742.79... 
## Restart number 15, logpost = -744.29... 
## Restart number 16, logpost = -745.03... 
## Restart number 17, logpost = -742.8... 
## Restart number 18, logpost = -745.04... 
## Restart number 19, logpost = -747.55... 
## Restart number 20, logpost = -745.03...

summary(fit3)

## __________________
## 
## Bayes-LCA
## Diagnostic Summary
## __________________
## 
## Hyper-Parameters: 
## 
##  Item Probabilities:
## 
##  alpha: 
##         Hallucination Activity Aggression Agitation Diurnal Affective
## Group 1             1        1          1         1       1         1
## Group 2             1        1          1         1       1         1
## Group 3             1        1          1         1       1         1
## 
##  beta: 
##         Hallucination Activity Aggression Agitation Diurnal Affective
## Group 1             1        1          1         1       1         1
## Group 2             1        1          1         1       1         1
## Group 3             1        1          1         1       1         1
## 
##  Class Probabilities:
## 
##  delta: 
## Group 1 Group 2 Group 3 
##       1       1       1 
## __________________
## 
## Method: EM algorithm  
## 
##  Number of iterations: 43 
## 
##  Log-Posterior Increase at Convergence: 0.0008143277 
## 
##  Log-Posterior: -742.7946 
## 
##  AIC: -1524.203 
## 
##  BIC: -1593.816

data(Alzheimer)
fit.gibbs <- blca.gibbs(Alzheimer, 2)

## Initialising sampler...starting burn-in.
## Burn-in completed...
## 1000 of 5000 samples completed...
## 2000 of 5000 samples completed...
## 3000 of 5000 samples completed...
## 4000 of 5000 samples completed...
## 5000 of 5000 samples completed...
## Sampling run completed.

## Warning in blca.gibbs(Alzheimer, 2): Label-switching (provisionally)
## corrected for - diagnostic plots are recommended. Use '?plot.blca' for
## details.

raftery.diag(as.mcmc(fit.gibbs))

## 
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
##                                                     
##               Burn-in  Total Lower bound  Dependence
##               (M)      (N)   (Nmin)       factor (I)
##  ClassProb 1  18       19341 3746         5.16      
##  ClassProb 2  32       29960 3746         8.00      
##  ItemProb 1 1 4        4955  3746         1.32      
##  ItemProb 1 2 8        10664 3746         2.85      
##  ItemProb 1 3 12       13382 3746         3.57      
##  ItemProb 1 4 18       19678 3746         5.25      
##  ItemProb 1 5 16       15640 3746         4.18      
##  ItemProb 1 6 16       19736 3746         5.27      
##  ItemProb 2 1 6        8244  3746         2.20      
##  ItemProb 2 2 15       18453 3746         4.93      
##  ItemProb 2 3 21       24489 3746         6.54      
##  ItemProb 2 4 18       20676 3746         5.52      
##  ItemProb 2 5 15       22893 3746         6.11      
##  ItemProb 2 6 21       27390 3746         7.31

#fit.gibbs <- blca.gibbs(Alzheimer, 2, K=100000, accept=0.1, burn.in=100)
plot(as.mcmc(fit.gibbs))

type1 <- c(0.8, 0.8, 0.2, 0.2)
type2 <- c(0.2, 0.2, 0.8, 0.8)
x<- rlca(1000, rbind(type1,type2), c(0.6,0.4))
set.seed(1)
fit <- blca(x, 2) ## EM algorithm used, warning returned

## Restart number 1, logpost = -2450.37... 
## New maximum found... Restart number 2, logpost = -2450.37... 
## Restart number 3, logpost = -2450.37... 
## New maximum found... Restart number 4, logpost = -2450.37... 
## Restart number 5, logpost = -2450.37...

print(fit) ## No posterior standard deviations returned

## 
## MAP Estimates:
##  
## 
## Item Probabilities:
##  
##         Col 1 Col 2 Col 3 Col 4
## Group 1  0.82 0.793 0.196 0.233
## Group 2  0.13 0.237 0.816 0.813
## 
## Membership Probabilities:
##  
## Group 1 Group 2 
##   0.595   0.405

## Warning: Posterior standard deviations not returned.

summary(fit)

## __________________
## 
## Bayes-LCA
## Diagnostic Summary
## __________________
## 
## Hyper-Parameters: 
## 
##  Item Probabilities:
## 
##  alpha: 
##         Col 1 Col 2 Col 3 Col 4
## Group 1     1     1     1     1
## Group 2     1     1     1     1
## 
##  beta: 
##         Col 1 Col 2 Col 3 Col 4
## Group 1     1     1     1     1
## Group 2     1     1     1     1
## 
##  Class Probabilities:
## 
##  delta: 
## Group 1 Group 2 
##       1       1 
## __________________
## 
## Method: EM algorithm  
## 
##  Number of iterations: 27 
## 
##  Log-Posterior Increase at Convergence: 0.001011169 
## 
##  Log-Posterior: -2450.371 
## 
##  AIC: -4918.742 
## 
##  BIC: -4962.912

set.seed(1)
fit2 <- blca(x, 2, method="em", sd=TRUE) ##No warning - same fit

## Restart number 1, logpost = -2450.37... 
## New maximum found... Restart number 2, logpost = -2450.37... 
## Restart number 3, logpost = -2450.37... 
## New maximum found... Restart number 4, logpost = -2450.37... 
## Restart number 5, logpost = -2450.37...

print(fit2) ##Posterior standard deviations returned

## 
## MAP Estimates:
##  
## 
## Item Probabilities:
##  
##         Col 1 Col 2 Col 3 Col 4
## Group 1  0.82 0.793 0.196 0.233
## Group 2  0.13 0.237 0.816 0.813
## 
## Membership Probabilities:
##  
## Group 1 Group 2 
##   0.595   0.405 
## 
## Posterior Standard Deviation Estimates:
##  
## 
## Item Probabilities:
##  
##         Col 1 Col 2 Col 3 Col 4
## Group 1 0.022 0.021 0.020 0.021
## Group 2 0.025 0.026 0.027 0.026
## 
## Membership Probabilities:
##  
## Group 1 Group 2 
##   0.023   0.023

set.seed(1)
##Variational Bayes approximation, with priors specified.
fit3 <- blca(x, 2, method="vb", delta=c(5,5), alpha=2, beta=1)

## Restart number 1, logpost = -6945.23...

print(fit3) ##Posterior standard deviations returned also.

## 
## MAP Estimates:
##  
## 
## Item Probabilities:
##  
##         Col 1 Col 2 Col 3 Col 4
## Group 1 0.821 0.794 0.197 0.234
## Group 2 0.135 0.240 0.815 0.813
## 
## Membership Probabilities:
##  
## Group 1 Group 2 
##   0.593   0.407 
## 
## Posterior Standard Deviation Estimates:
##  
## 
## Item Probabilities:
##  
##         Col 1 Col 2 Col 3 Col 4
## Group 1 0.016 0.017 0.016 0.017
## Group 2 0.017 0.021 0.019 0.019
## 
## Membership Probabilities:
##  
## Group 1 Group 2 
##   0.015   0.015

par(mfrow=c(3,2))
plot(fit3, which=3:4)
par(mfrow=c(1,1))