In this example, we write a function to augment the default output of randomLCA() to make it easier for the analyst to interpret the results.
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 will write a function wrapper for randomLCA that does some additional work in a generic fashion. This will allow easier estimation of other models. We annotate the function to explain what we're doing. The resulting objects are outcomep, which contains the outcome probabilities, and classp, with the class probabilities.
runlca = function(df, nclass=2, names=c("item"), verbose=FALSE) {
nvars = dim(df)[2]
# create a list of names for the items
if (length(names)==1) { names = rep(names, nvars) }
# include only complete cases
bigtable = table(na.omit(df))
allpatterns = as.data.frame(ftable(bigtable))
# keep only the patterns that occur
nonzeropatterns = allpatterns[allpatterns$Freq > 0,]
# fit the model
results = randomLCA(nonzeropatterns[,1:nvars],
nonzeropatterns$Freq, nclass=nclass, calcSE=FALSE)
# display available sample size
cat("nobs=", results$nobs, "\n")
oldopt = options(digits=2)
if (verbose==TRUE) { # display patterns
whichclass = apply(results$classprob, 1, which.max)
nonzeropatterns$class = whichclass
print(nonzeropatterns[order(whichclass),])
}
print(summary(results))
resvals = cbind(results$outcomep, results$classp)
# label the margins with our desired variable names
# (plus class probability)
colnames(resvals) = c(names, "classprob")
# annotate standard output with rounded values
print(round(resvals, 2))
options(oldopt)
return(results)
}
Now let's apply the function. 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 to the function. Then we call the runlca() function. By specifying the verbose option the code displays each of the patterns, sorted by which class it is in (based on the highest predicted probability).
cesdcut = ifelse(cesd>20, 1, 0)
smallds = data.frame(homeless, cesdcut, satreat, linkstatus)
results = runlca(smallds, nclass=3,
names=c("homeless", "cesd", "satreat", "linkstatus"),
verbose=TRUE)
This generates the following output:
nobs= 431
homeless cesdcut satreat linkstatus Freq class
5 0 0 1 0 16 1
7 0 1 1 0 33 1
6 1 0 1 0 4 2
8 1 1 1 0 37 2
13 0 0 1 1 1 2
14 1 0 1 1 4 2
15 0 1 1 1 9 2
16 1 1 1 1 23 2
1 0 0 0 0 17 3
2 1 0 0 0 15 3
3 0 1 0 0 82 3
4 1 1 0 0 64 3
9 0 0 0 1 10 3
10 1 0 0 1 9 3
11 0 1 0 1 62 3
12 1 1 0 1 45 3
Classes AIC BIC logLik
3 2093 2150 -1032
Class probabilities
Class 1 Class 2 Class 3
0.07846 0.21621 0.70534
Outcome probabilities
homeless cesd satreat linkstatus classprob
[1,] 0.00 0.58 1 0.00 0.08
[2,] 0.73 0.88 1 0.40 0.22
[3,] 0.44 0.83 0 0.41 0.71
The results are equivalent to the results from the prior example, but the predicted classes are listed, and the class probabilities (and proportion endorsing the item) are more clearly discernible. It might be useful in a later iteration of the function to add some blank lines and the proportion of the seeds that resulted in the maximum likelihood.
1 comment:
I'm the author of randomLCA, thanks for your comments. I've cleaned up the output a fair bit and made some other changes. They should be available on CRAN now.
Post a Comment