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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
|
`LRtest.Rm` <-
function(object, splitcr = "median", se = TRUE)
{
# performs Andersen LR-test
# object... object of class RM
# splitcr... splitting criterion for LR-groups. "all.r" corresponds to a complete
# raw score split (r=1,...,k-1), "median" to a median raw score split,
# "mean" corresponds to the mean raw score split.
# optionally also a vector of length n for group split can be submitted.
# se...whether standard errors should be computed
call<-match.call()
spl.gr<-NULL
X.original<-object$X
if((length(splitcr) > 1) & is.character(splitcr)){ # if splitcr is character vector, treated as factor
splitcr<-as.factor(splitcr)
}
if(is.factor(splitcr)){
spl.nam<-deparse(substitute(splitcr))
spl.lev<-levels(splitcr)
spl.gr<-paste(spl.nam,spl.lev,sep=" ")
splitcr<-unclass(splitcr)
}
numsplit<-is.numeric(splitcr)
if (any(is.na(object$X))) {
if (!numsplit && splitcr=="mean") { #mean split
spl.gr<-c("Raw Scores < Mean", "Raw Scores >= Mean")
X<-object$X
# calculates index for NA groups
# from person.parameter.eRm
dichX <- ifelse(is.na(X),1,0)
strdata <- apply(dichX,1,function(x) {paste(x,collapse="")})
gmemb <- as.vector(data.matrix(data.frame(strdata)))
gindx<-unique(gmemb)
rsum.all<-rowSums(X,na.rm=T)
grmeans<-tapply(rsum.all,gmemb,mean) #sorted
ngr<-table(gmemb) #sorted
m.all<-rep(grmeans,ngr) #sorted,expanded
rsum.all<-rsum.all[order(gmemb)]
spl<-ifelse(rsum.all<m.all,1,2)
splitcr<-spl
object$X<-X[order(gmemb),]
}
if (!numsplit && splitcr=="median") { #median split
spl.gr<-c("Raw Scores <= Median", "Raw Scores > Median")
# cat("Warning message: Persons with median raw scores are assigned to the lower raw score group!\n")
X<-object$X
# calculates index for NA groups
# from person.parameter.eRm
dichX <- ifelse(is.na(X),1,0)
strdata <- apply(dichX,1,function(x) {paste(x,collapse="")})
gmemb <- as.vector(data.matrix(data.frame(strdata)))
gindx<-unique(gmemb)
rsum.all<-rowSums(X,na.rm=T)
grmed<-tapply(rsum.all,gmemb,median) #sorted
ngr<-table(gmemb) #sorted
m.all<-rep(grmed,ngr) #sorted,expanded
rsum.all<-rsum.all[order(gmemb)]
spl<-ifelse(rsum.all<=m.all,1,2)
splitcr<-spl
object$X<-X[order(gmemb),]
}
}
if (!is.numeric(splitcr)) {
if (splitcr=="all.r") { #full raw score split ### begin MjM 2012-03-18
rvind <- rowSums(object$X, na.rm=TRUE) #person raw scoobject
excl_0_k <- (rvind > 0) & (rvind < sum(apply(object$X, 2, max, na.rm=T)))
Xlist <- by(object$X[excl_0_k,], rvind[excl_0_k], function(x) x)
names(Xlist) <- as.list(paste("Raw Score =", sort(unique(rvind[excl_0_k]))))
spl.gr <- unlist(names(Xlist))
} ### end MjM 2012-03-18
if (splitcr=="median") { #median split
spl.gr<-c("Raw Scores <= Median", "Raw Scores > Median")
#removed rh 2010-12-17
#cat("Warning message: Persons with median raw scores are assigned to the lower raw score group!\n")
rv <- apply(object$X,1,sum,na.rm=TRUE)
rvsplit <- median(rv)
rvind <- rep(0,length(rv))
rvind[rv > rvsplit] <- 1 #group with highraw scoobject
Xlist <- by(object$X,rvind,function(x) x)
names(Xlist) <- list("low","high")
}
if (splitcr=="mean") { #mean split
spl.gr<-c("Raw Scores < Mean", "Raw Scores >= Mean")
rv <- apply(object$X,1,sum,na.rm=TRUE)
rvsplit <- mean(rv)
rvind <- rep(0,length(rv))
rvind[rv > rvsplit] <- 1 #group with highraw scoobject
Xlist <- by(object$X,rvind,function(x) x)
names(Xlist) <- list("low","high")
}
}
if (is.numeric(splitcr)) { #manual raw score split
spl.nam<-deparse(substitute(splitcr))
if (length(splitcr)!=dim(object$X)[1]) stop("Mismatch between length of split vector and number of persons!")
if (any(is.na(splitcr))) stop("Split vector should not contain NA's")
rvind <- splitcr
Xlist <- by(object$X,rvind, function(x) x)
names(Xlist) <- as.list(sort(unique(splitcr)))
if(is.null(spl.gr)){
spl.lev<-names(Xlist)
spl.gr<-paste(spl.nam,spl.lev,sep=" ")
}
}
#----------item to be deleted---------------
del.pos.l <- lapply(Xlist, function(x) {
it.sub <- datcheck.LRtest(x,object$X,object$model) #items to be removed within subgroup
})
del.pos <- unique(unlist(del.pos.l))
if (length(del.pos) >= (ncol(object$X)-1)) {
stop("\nNo items with appropriate response patterns left to perform LR-test!\n")
}
if(length(del.pos) > 0){ ### begin MjM 2013-01-27
warning(paste0(
"\n",
prettyPaste("The following items were excluded due to inappropriate response patterns within subgroups:"),
"\n",
paste(colnames(object$X)[del.pos], collapse=" "),
"\n\n",
prettyPaste("Full and subgroup models are estimated without these items!")
), immediate.=TRUE)
} ### end MjM 2013-01-27
if (length(del.pos) > 0) {
X.el <- object$X[,-(del.pos)]
} else {
X.el <- object$X
}
if(ifelse(length(splitcr) == 1, splitcr != "all.r", TRUE)){ ### begin MjM 2012-03-18 # for all cases except "all.r"
Xlist.n <- by(X.el, rvind, function(y) y)
names(Xlist.n) <- names(Xlist)
if (length(del.pos) > 0) Xlist.n <- c(Xlist.n,list(X.el)) # X.el added since we must refit whole group without del.pos items
} else {
Xlist.n <- by(X.el[excl_0_k,], rvind[excl_0_k], function(y) y)
names(Xlist.n) <- names(Xlist)
Xlist.n <- c(Xlist.n,list(X.el[excl_0_k,])) # X.el added since we must refit whole group without del.pos items
} ### end MjM 2012-03-18
if (object$model=="RM") {
likpar <- sapply(Xlist.n,function(x) { #matrix with loglik and npar for each subgroup
objectg <- RM(x,se=se)
likg <- objectg$loglik
nparg <- length(objectg$etapar)
# betalab <- colnames(objectg$X)
list(likg,nparg,objectg$betapar,objectg$etapar,objectg$se.beta,outobj=objectg) # rh outobj added
###list(likg,nparg,objectg$betapar,objectg$etapar,objectg$se.beta) # rh outobj added
})
}
if (object$model=="PCM") {
likpar <- sapply(Xlist.n,function(x) { #matrix with loglik and npar for each subgroup
objectg <- PCM(x,se=se)
likg <- objectg$loglik
nparg <- length(objectg$etapar)
list(likg,nparg,objectg$betapar,objectg$etapar,objectg$se.beta,outobj=objectg) # rh outobj added
###list(likg,nparg,objectg$betapar,objectg$etapar,objectg$se.beta) # rh outobj added
})
}
if (object$model=="RSM") {
likpar <- sapply(Xlist.n,function(x) { #matrix with loglik and npar for each subgroup
objectg <- RSM(x,se=se)
likg <- objectg$loglik
nparg <- length(objectg$etapar)
list(likg,nparg,objectg$betapar,objectg$etapar,objectg$se.beta,outobj=objectg) # rh outobj added
###list(likg,nparg,objectg$betapar,objectg$etapar,objectg$se.beta) # rh outobj added
})
}
## extract fitted splitgroup models # rh 02-05-2010
if(ifelse(length(splitcr) == 1, splitcr != "all.r", TRUE)){ ### begin MjM 2012-03-18
fitobj <- likpar[6, 1:length(unique(rvind))]
} else {
fitobj <- likpar[6, 1:length(unique(rvind[excl_0_k]))]
} ### end MjM 2012-03-18
likpar <- likpar[-6,]
if((length(del.pos) > 0) | ifelse(length(splitcr) == 1, splitcr == "all.r", FALSE)) { #re-estimate full model ### MjM 2012-03-18
pos <- length(Xlist.n) #position of the full model
loglik.all <- likpar[1,pos][[1]] #loglik full model
# etapar.all <- rep(0,likpar[2,pos]) #etapar full model (filled with 0 for df computation)
etapar.all <- rep(0, unlist(likpar[2,pos])) #etapar full model (filled with 0 for df computation)
likpar <- likpar[,-pos]
Xlist.n <- Xlist.n[-pos]
} else {
loglik.all <- object$loglik
etapar.all <- object$etapar
}
loglikg <- sum(unlist(likpar[1,])) #sum of likelihood value for subgroups
LR <- 2*(abs(loglikg-loglik.all)) #LR value
df = sum(unlist(likpar[2,]))-(length(etapar.all)) #final degrees of freedom
pvalue <- 1 - pchisq(LR, df) #pvalue
betalist <- likpar[3,] #organizing betalist
result <- list(X=X.original, X.list=Xlist.n, model=object$model,LR=LR,
df=df, pvalue=pvalue, likgroup=unlist(likpar[1,],use.names=FALSE),
betalist=betalist, etalist=likpar[4,],selist=likpar[5,], spl.gr=spl.gr, call=call, fitobj=fitobj) ## rh fitobj added
class(result) <- "LR"
return(result)
}
|
