1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
##' Randomly sample response patterns given a list of items
##'
##' Returns a random sample of response patterns given a list of item
##' models and parameters. If \code{grp} is given then theta, items, params,
##' mean, and cov can be omitted.
##'
##' @template detail-group
##' @name rpf.sample
##' @param theta either a vector (for 1 dimension) or a matrix (for >1
##' dimension) of person abilities or the number of response patterns
##' to generate randomly
##' @param items a list of item models
##' @param params a list or matrix of item parameters. If omitted, random item
##' parameters are generated for each item model.
##' @template arg-dots
##' @param prefix Column names are taken from param or items.
##' If no column names are available, some will be generated using
##' the given prefix.
##' @param mean mean vector of latent distribution (optional)
##' @param cov covariance matrix of latent distribution (optional)
##' @param mcar proportion of generated data to set to NA (missing completely at random)
##' @template arg-grp
##' @return Returns a data frame of response patterns
##' @examples
##' # 1 dimensional items
##' i1 <- rpf.drm()
##' i1.p <- rpf.rparam(i1)
##' i2 <- rpf.nrm(outcomes=3)
##' i2.p <- rpf.rparam(i2)
##' rpf.sample(5, list(i1,i2), list(i1.p, i2.p))
##' @seealso \code{\link{sample}}
rpf.sample <- function(theta, items, params, ..., prefix="i",
mean=NULL, cov=NULL, mcar=0.0, grp=NULL)
{
garbageArguments <- list(...)
if (length(garbageArguments) > 0) {
stop("rpf.sample does not accept values for the '...' argument")
}
if (!missing(grp)) {
if (missing(items)) items <- grp$spec
if (missing(params) && !is.null(grp$param)) params <- grp$param
if (missing(mean) && !is.null(grp$mean)) mean <- grp$mean
if (missing(cov) && !is.null(grp$cov)) cov <- grp$cov
if (missing(theta) && !is.null(grp$data)) theta <- nrow(grp$data)
}
numItems <- length(items)
maxDim <- max(vapply(items, function(i) i@factors, 0))
if (maxDim > 1) {
design <- matrix(rep(1:maxDim, numItems), nrow=maxDim)
design[sapply(items, function(i) 1:maxDim > i@factors)] <- NA
} else {
design <- matrix(rep(1, numItems), nrow=1)
}
maxAbilities <- max(design, na.rm=TRUE)
numPeople <- NA
if (is.numeric(theta) && length(theta) == 1) {
if (theta <= 1) stop("Request at least 2 samples")
numPeople <- theta
if (missing(mean)) mean <- rep(0, maxAbilities)
if (length(mean) != maxAbilities) stop(paste("Mean vector must have length",maxAbilities,"not",length(mean)))
if (missing(cov)) cov <- diag(maxAbilities)
if (any(dim(cov) != maxAbilities)) stop(paste("Cov matrix must be square matrices of size",maxAbilities))
theta <- array(t(mvtnorm::rmvnorm(numPeople, mean=mean, sigma=cov)),
dim=c(maxAbilities, numPeople))
} else if (maxDim == 1 && is.vector(theta)) {
numPeople <- length(theta)
theta <- array(theta, dim=c(maxAbilities, numPeople))
} else {
if (dim(theta)[1] > maxAbilities) {
stop(paste("Only", maxAbilities, "abilities but theta provides", dim(theta)[1],
"-- maybe transpose theta?"))
}
numPeople <- dim(theta)[2]
}
if (missing(params)) {
params <- lapply(items, rpf.rparam)
}
outcomes <- vapply(items, function(i) i@outcomes, 0)
name <- colnames(params)
if (is.null(name)) name <- names(items)
maxParam <- max(sapply(items, rpf.numParam))
if (is.matrix(params) && nrow(params) > maxParam) {
warning(paste("Item parameter matrix has", nrow(params) - maxParam, "extra rows"))
}
ret <- list()
for (ix in 1:numItems) {
i <- items[[ix]]
param <- c()
if (is.list(params)) {
param <- params[[ix]]
} else {
param <- params[,ix] # item parameters are in columns
}
if (length(param) < rpf.numParam(i)) {
stop(paste("Item",class(i),"needs",rpf.numParam(i),
"parameters but only",length(param),"given"))
}
cols <- design[,ix]
cols <- cols[!is.na(cols)]
i.theta <- theta[cols,,drop=FALSE]
P <- rpf.prob(i, param[1:rpf.numParam(i)], i.theta)
# if (any(is.na(P))) stop(paste("Item", i@spec, "with param", param," produced NAs"))
ret1 <- apply(P, 2, sample, x=1:i@outcomes, size=1, replace=F)
labels <- 1:i@outcomes
if (!missing(grp) && !is.null(grp$data)) labels <- levels(grp$data[1, name[ix]])
ret1 <- factor(ret1, levels=1:i@outcomes, ordered=TRUE, labels=labels)
attr(ret1, 'mxFactor') <- TRUE # for OpenMx
ret[[ix]] <- ret1
}
ret <- as.data.frame(ret, optional=TRUE)
if (is.null(name)) name <- paste(prefix,1:numItems,sep="")
colnames(ret) <- name
if (mcar > 0) {
size <- prod(dim(ret))
mask <- rep(FALSE, size)
mask[sample.int(size, size * mcar)] <- TRUE
shaped.mask <- array(mask, dim=dim(ret))
ret[shaped.mask] <- NA
}
return(ret)
}
|
