In this example, we fit the same model using the randomLCA() function within the package of the same name.
R
We begin by reading in the data.
ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
attach(ds)
library(randomLCA)
We start by creating a dichotomous variable with high scores on the CESD, and put this together as part of a dataframe to be given as input.
cesdcut = ifelse(cesd>20, 1, 0)
smallds = na.omit(data.frame(homeless, cesdcut, satreat, linkstatus))
results = randomLCA(smallds, nclass=3, notrials=1000)
summary(results)
This generates the following output:
Classes AIC BIC logLik
3 2092.968 2149.893 -1032.484
Class probabilities
Class 1 Class 2 Class 3
0.07846 0.70534 0.21620
Outcome probabilities
homeless cesdcut satreat linkstatus
Class 1 9.465e-06 0.5786 1.000e+00 9.538e-06
Class 2 4.375e-01 0.8322 9.988e-06 4.145e-01
Class 3 7.297e-01 0.8846 1.000e+00 3.971e-01
The results are equivalent to the results from the prior example, though the scientific notation for the observed prevalences in each class are hard to read. Other objects are available from the returned value, though they are also not in an easily digestible form:
> names(results)
[1] "fit" "nclass" "classp" "outcomep"
[5] "se" "np" "nobs" "logLik"
[9] "observed" "fitted" "deviance" "classprob"
[13] "bics" "random" "level2" "byclass"
[17] "blocksize" "call" "probit" "quadpoints"
[21] "patterns" "notrials" "freq"
> results$patterns
homeless cesdcut satreat linkstatus
1 0 0 0 0
2 0 0 0 1
3 0 0 1 0
4 0 0 1 1
5 0 1 0 0
6 0 1 0 1
7 0 1 1 0
8 0 1 1 1
9 1 0 0 0
10 1 0 0 1
11 1 0 1 0
12 1 0 1 1
13 1 1 0 0
14 1 1 0 1
15 1 1 1 0
16 1 1 1 1
> results$freq
[1] 17 10 16 1 82 62 33 9 15 9 4 4 64 45 37 23
We'll address workarounds for these shortcomings in a future entry.
No comments:
Post a Comment