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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
plotDIF <- function(object, item.subset=NULL, gamma = 0.95, main=NULL,
xlim=NULL, xlab=" ", ylab=" ", col=NULL, distance,
splitnames=NULL, leg=FALSE, legpos="bottomleft", ...){
if("LR" %in% class(object)){
object <- list(object)
} else if(is.list(object)) {
checklr <- sapply(object,class)
if(!all(checklr=="LR")) stop("Elements of '",deparse(substitute(object)), "' must must be LRtest objects!")
} else if(!is.list(object)) {
stop(deparse(substitute(object)), "must be a list of LRtest objects!")
}
# extract number of LRtest objects
M <- length(sapply(object, function(x) length(x)))
# Confidence plot only for LRtest objects
for(p in 1:M){
if((object[[p]]$model == "LLTM") || (object[[p]]$model == "LRSM") || (object[[p]]$model == "LPCM")){
stop("Confidence Plot is computed only for LRtest objects (RM, PCM, RSM)!")
} else if(is.na(sum(unlist(object[[p]]$selist))) == TRUE){
stop("Confidence Plot is computed only for LRtest objects (RM) with standard errors (se=TRUE)!")
}
}
# confidences list for storing confints
confidences1 <- vector("list")
# for labeling list entries
nam.lab <- vector("character")
# subgroups splits
n2 <- sapply(object, function(x) length(x$spl.gr))
# loops for computing thresholds on LRtest objects
for(m in 1:M){ # confidences for dichotomous items
if(object[[m]]$model == "RM"){
confidences1[[m]] <- lapply(object[[m]]$fitobj,function(x){-confint(x, level=gamma)})
} else { # confidences for polytomous items
confidences1[[m]] <- lapply(object[[m]]$fitobj, function(x){confint(thresholds(x),level=gamma)})
}
}
if(is.null(names(object)) == TRUE){
names(confidences1) <- paste("LRtest", 1:M, sep="")
} else {
names(confidences1) <- names(object)
}
confidences <- do.call(c,lapply(confidences1,function(x) x[1:length(x)]))
if(missing(distance)) distance <- .7/(length(confidences))
if((distance <= 0) | (distance >= .5)) stop("distance must not be >= .5 or <= 0")
model.vec <- vector("character")
for(p in 1:M) model.vec[p] <- object[[p]]$model
if(any(model.vec == "PCM") || any(model.vec == "RSM")){
model <- "PCM"
} else {
model <- "RM"
}
# extracting the longest element of confidences for definition of tickpositions and ticklabels
# (snatches at the confidences-object index)
factorlist <- (unique(unlist(lapply(confidences, function(x) dimnames(x)[[1]]))))
#maxlist <- max(order(factorlist))
if(is.null(item.subset)){
if(model == "PCM"){
y.lab <- sub("thresh beta ", "", factorlist)
} else {
y.lab <- sub("beta ", "", factorlist)
}
} else if(is.character(item.subset)){ # item subset specified as character
if(model == "PCM"){
y.lab <- sub("(.+)[.][^.]+$", "\\1", sub("thresh beta ", "", factorlist)) # search only for a "." separation after item label
categ <- gsub("^.*\\.(.*)$","\\1", factorlist) # extract item categories - search only for a "." separation
y.lab.id <- y.lab %in% item.subset
y.lab1 <- y.lab[y.lab.id]
categ1 <- categ[y.lab.id]
y.lab <- paste(y.lab1, categ1, sep=".") # stick item categories and names together again
factorlist <- factorlist[y.lab.id]
} else {
y.lab <- sub("beta ", "", factorlist) # search only for a "." separation after item label
y.lab.id <- y.lab %in% item.subset
y.lab <- y.lab[y.lab.id]
factorlist <- factorlist[y.lab.id]
}
} else { # item subset specified as position number (index in data matrix)
if(model == "PCM"){
y.lab <- sub("(.+)[.][^.]+$", "\\1", sub("thresh beta ", "", factorlist)) # search only for a "." separation after item label
categ <- gsub("^.*\\.(.*)$","\\1", factorlist) # extract item categories - search only for a "." separation
y.lab2 <- unique(y.lab)[item.subset]
y.lab.id <- y.lab %in% y.lab2
y.lab1 <- y.lab[y.lab.id]
categ1 <- categ[y.lab.id]
y.lab <- paste(y.lab1, categ1, sep=".") # stick item categories and names together again
factorlist <- factorlist[y.lab.id]
} else {
y.lab <- sub("beta ", "", factorlist)
y.lab2 <- unique(y.lab)[item.subset]
y.lab.id <- y.lab %in% y.lab2
y.lab <- y.lab[y.lab.id]
factorlist <- factorlist[y.lab.id]
}
}
# setting range of xaxis
if(is.null(xlim)){ xlim <- range(unlist(confidences)) }
# setting tickpositions
tickpos <- 1:(length(y.lab)) # + 3.5/distance mm 2011-06-03
lty <- unlist(lapply(n2, function(i)1:i))
# defining the plot
if(is.null(main)){ main<-paste("Confidence plot") }
plot(xlim, xlim=xlim, ylim=c(1,length(factorlist))+c(-.5,+.5), type="n", yaxt="n", main=main, xlab=xlab,ylab=ylab,...) # rh 2011-03-23 reverse ylim added
axis(2, at=tickpos, labels=y.lab, cex.axis=0.7, las=2)
if(is.null(col)){
for(k in 1:length(confidences)) {
for(l in 1:length(factorlist)) {
lines(as.data.frame(confidences[[k]])[factorlist[l],],
rep(seq(l-.5+distance, l+.5-distance, length.out=length(confidences))[k], 2),
type="b", pch=20, col=length(cumsum(n2)[cumsum(n2) < k])+1, lty=lty[k])
}
}
} else {
col <- rep(col, n2)
for(k in 1:length(confidences)) {
for(l in 1:length(factorlist)) {
lines(as.data.frame(confidences[[k]])[factorlist[l],],
rep(seq(l-.5+distance, l+.5-distance, length.out=length(confidences))[k], 2),
type="b", pch=20, col=col[k], lty=lty[k])
}
}
}
# doing nicer legend labels
if(is.null(splitnames)==FALSE){ names(confidences) <- splitnames }
if(leg == TRUE){
linespread <- .7 + .3 * (1/length(confidences))
if(is.null(col)){ #col <- rep(1:length(n2), each=n2) rh 2011-03-18
col <- rep(1:length(n2), n2) # legend(legpos, rev(paste(names(confidences))), col=rev(col), lty=rev(lty))
legend(legpos, rev(paste(names(confidences))), y.intersp=linespread, col=rev(col), lty=rev(lty))
} else { # legend(legpos, rev(paste(names(confidences))), col=rev(col), lty=rev(lty))
legend(legpos, rev(paste(names(confidences))), y.intersp=linespread, col=rev(col), lty=rev(lty))
}
}
invisible(list(confints=confidences1)) #rh 2011-03-18
}
|
