Nothing
## File Name: tam_A_matrix2.R
## File Version: 9.171
## function tam_A_matrix2
tam_A_matrix2__INACTIVE <- function( resp, formulaA=~ item + item*step, facets=NULL,
constraint=c("cases", "items"), progress=FALSE, maxKi=NULL, Q=Q )
{
z0 <- Sys.time()
### redefine facets matrix
facets0 <- facets
NF <- length(facets)
facet.list <- as.list( 1:NF )
names(facet.list) <- colnames(facets)
if (NF==0){ facet.list <- NULL }
if (NF>0){
for (ff in 1:NF){
# ff <- 2
uff <- sort( unique( facets[,ff] ) )
facets[,ff] <- match( facets[,ff], uff )
facet.list[[ff]] <- data.frame(
"facet.label"=paste0( colnames(facets)[ff], uff ),
"facet.index"=paste0( colnames(facets)[ff], seq(1,length(uff) ) ) )
}
}
#cat(" +++ v62" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
### Basic Information and Initializations
constraint <- match.arg(constraint)
if ( is.null(maxKi) ){
maxKi <- apply( resp, 2, max, na.rm=TRUE )
}
maxK <- max( maxKi )
nI <- ncol(resp)
# cat(" +++ v70" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
# stop processing if there are items with a maximum score of 0
i11 <- names(maxKi)[ maxKi==0 ]
if ( length(i11) > 0 ){
stop( cat( "Items with maximum score of 0:", paste(i11, collapse=" " ) ) )
}
tf <- stats::terms( formulaA )
fvars <- as.vector( attr(tf,"variables"), mode="character" )[-1]
#cat("fvars 212") ; print(fvars)
otherFacets <- setdiff( fvars, c("item", "step") )
contr.list <- as.list( rep( "contr.sum", length(fvars) ) )
names( contr.list ) <- fvars
#cat(" +++ v80" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
#****
# ARb: 2013-03-27
# no contrasts for items
nitems <- ncol(resp)
# contr.list[["item"]] <- diag(1,nitems)
# TK: 2014-03-12
# Consider long vector response matrix and "item" in facets input
expand.list <- vector(mode="list", length=0)
if( "item" %in% fvars ) expand.list <- c(expand.list, if("item" %in% names(facet.list)) list(as.factor(sort(unique(facets[,"item"])))) else list(factor(1:nI)) )
if( "step" %in% fvars ) expand.list <- c(expand.list, if("step" %in% names(facet.list)) list(as.factor(sort(unique(facets[,"step"])))) else list(factor(1:maxK)) )
if( length( otherFacets )==1) expand.list <- c(expand.list, list(factor(1:max(facets[, otherFacets]))) )
if( length( otherFacets ) > 1 ) expand.list <- c(expand.list, sapply( otherFacets, FUN=function(ff) as.factor(1:max(facets[, ff])), simplify=FALSE ))
names( expand.list ) <- fvars
#cat(" +++ v100" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
# expand.list <- expand.list[ !unlist( lapply(expand.list, is.null) ) ]
for (vv in seq(1, length(expand.list) ) ){
expand.list[[vv]] <- paste( expand.list[[vv]] )
}
# cat(" +++ v110" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
g2 <- g1 <- expand.grid(expand.list)
diffK <- ( stats::sd( maxKi) > 1e-10 )
# diffK <- FALSE
diffK <- TRUE
# reduced combinations of items
if (diffK){
I <- length(maxKi)
# g1 <- as.data.frame(g1)
for (ii in 1:I){
ind <- which( ( ( as.numeric(paste0(g1$item))==ii ) &
( as.numeric(paste0(g1$step)) > maxKi[ii] ) ) )
if ( length(ind) > 0 ){
g1 <- g1[ - ind, ]
}
}
}
X <- rownames.design2( g1 )
if (diffK){
X2 <- rownames.design2( g2 )
}
#**** include maxKi (2014-05-30)
# X$maxKi <- maxKi[ X$item ]
# X$par.exist <- 1 * ( as.numeric(paste(X$step)) <=X$maxKi )
#****
### constraints and formulaA
if( constraint=="cases" ) formulaA <- stats::update.formula(formulaA, ~0+.)
NX <- ncol(X)
for (ff in 1:NX){
uff <- length( unique(X[,ff] ) )
if (uff==1){ cat(paste0(" - facet ",
colnames(X)[ff], " does only have one level!" ), "\n") }
}
# cat(" +++ v120" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
mm <- - stats::model.matrix(formulaA, X, contrasts=contr.list)
if (diffK){
mm2 <- - stats::model.matrix(formulaA, X2, contrasts=contr.list)
}
# mm <- - model.matrix(formulaA, X )
if( constraint=="items" ){ mm <- mm[,-1] }
############################################################
###*** ARb 2013-03-28
### generate all interactions
xsi.constr <- .generate.interactions2(X, facets, formulaA, mm )
###############################################################
#cat(" +++ v130" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
### Postprocessing
# model.matrix _ case: step in fvars
if( "step" %in% fvars ){
if( ncol( attr(tf, "factors") )==1 ){
return( warning("Can't proceed the estimation:
Factor of order 1 other than step must be specified.") )
}
if( all( attr(tf, "factors")["step",] !=1 ) ){
return( warning("Can't proceed the estimation:
Lower-order term is missing.") )
}
A <- NULL
stepgroups <- unique( gsub( "(^|-)+step([[:digit:]])*", "\\1step([[:digit:]])*",
x=rownames(X), perl=TRUE ) )
X.out <- data.frame(as.matrix(X), stringsAsFactors=FALSE)
#cat(" +++ v150" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
if (progress){
cat(" o Create design matrix A\n")
ip <- length(stepgroups)
VP <- min( ip, 10 )
cat(paste0(" |",paste0( rep("*", VP), collapse=""), "|\n"))
cat(" |") ; utils::flush.console()
if (VP<10){ disp_progress <- 1:ip } else {
disp_progress <- 100* ( 1:ip ) / (ip+1)
disp_progress <- sapply( seq(5,95,10), FUN=function(pp){ # pp <- 5
which.min( abs( disp_progress - pp ) )[1] }
)
}
}
#******
# collect xsi parameters to be excluded
xsi.elim.index <- xsi.elim <- NULL
ii <- 0 ; vv <- 1
for( sg in stepgroups ){
mm1 <- grep(paste0("(", sg, ")+$"), rownames(mm))
ind2 <- grep( paste0("(", sg, ")+$"), rownames(mm))
mm.sg.temp <- rbind( 0, apply( mm[ ind2,,drop=FALSE], 2, cumsum ) )
if ( is.null(rownames(mm.sg.temp)) ){
rownames(mm.sg.temp) <- paste0("rn", seq(0,nrow(mm.sg.temp)-1) )
}
# substitute the following line later if ...
rownames(mm.sg.temp)[1] <- gsub("step([[:digit:]])*", "step0", sg, fixed=TRUE)
# rownames(mm.sg.temp)[1] <- gsub("step([[:digit:]])*", "step0", sg, fixed=FALSE, perl=TRUE)
rownames(mm.sg.temp)[-1] <- rownames(mm[ind2,,drop=FALSE])
#****
# set entries to zero if there are no categories in data
sg1 <- strsplit( sg, split="-")[[1]]
ii <- as.numeric( gsub("item", "", sg1[1] ) )
maxKi_ii <- maxKi[ii]
if (is.na(ii)){
maxKi_ii <- maxK
}
if ( maxKi_ii < maxK ){
for (kk in (maxKi[ii]+1):maxK){
# kk <- 2
# set rows in A matrix to zero
mm.sg.temp[ grep( paste0( "-step", kk ), rownames(mm.sg.temp) ), ] <- NA
# mm.sg.temp[ grep( paste0( "-step", kk ), rownames(mm.sg.temp) ), ] <- 0
i1 <- grep( paste0(sg1[1],"\\:" ), colnames(mm.sg.temp), value=TRUE)
i2 <- grep( paste0(":step", kk-1), colnames(mm.sg.temp), value=TRUE)
i3 <- intersect( i1, i2 )
if ( length(i3) > 0 ){
xsi.elim <- c( xsi.elim, i3 )
xsi.elim.index <- c( xsi.elim.index,
which( colnames(mm.sg.temp ) %in% i3 ) )
# mm.sg.temp[, i3 ] <- NA
# mm.sg.temp[ ! ( is.na( mm.sg.temp[, i3 ] ) ), i3 ] <- 0
#@@@ suggestion Michal Modzelewski
mm.sg.temp[, i3 ][ ! ( is.na( mm.sg.temp[, i3 ] ) )] <- 0
} #**** end i3
#*******************
}
}
A <- rbind(A, mm.sg.temp)
if ( maxKi_ii < maxK ){
for (kk in (maxKi[ii]+1):maxK){
vv <- paste0( sg1[1], "-step", kk )
X.out[ grep( vv, rownames(X.out) ), 2 ] <- 0
}
}
x.sg.temp <- X.out[grep(sg, rownames(X.out))[1], ]
x.sg.temp[,"step"] <- 0
rownames(x.sg.temp) <- gsub("step([[:digit:]])*", "step0", sg, fixed=TRUE)
X.out <- rbind(X.out, x.sg.temp)
progress <- progress & ( ! is.na(ii) )
if ( progress ){
ii <- ii+1
if (( ii==disp_progress[vv] ) & (vv<=10) ){
cat("-") ; utils::flush.console()
vv <- vv+1
}
} # end progress
} # end stepgroups
if ( progress ){
cat("|\n") ; utils::flush.console()
}
# cat(" +++ v160" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
} else {
# model.matrix _ case: step not in fvars
rownames(mm) <- paste( rownames(X), "-step1", sep="")
A <- mm
for( kk in setdiff(0:maxK, 1) ){
mm.k.temp <- mm*kk
rownames(mm.k.temp) <- paste( rownames(X), "-step", kk, sep="")
A <- rbind(A, mm.k.temp)
}
#cat(" +++ v170" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
X.out <- expand.grid( c( expand.list, list("step"=factor(0:maxK)) ) )
X.out <- rownames.design2( data.frame(as.matrix(X.out), stringsAsFactors=FALSE) )
}# end step in fvars
#***
# set entries in A to zero for constraints
A[ rowMeans(is.na(A)) < 1, xsi.elim ] <- 0
#cat(" +++ v180" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
# facet design
facet.design <- list( "facets"=facets, "facets.orig"=facets0,
"facet.list"=facet.list[otherFacets])
A <- A[ ! duplicated( rownames(A) ), ]
if ( max(apply(resp,2,max,na.rm=TRUE)) > 9 ){
A <- A[order(rownames(A)),,drop=FALSE]
}
X.out <- X.out[order(rownames(X.out)),,drop=FALSE]
#*** elimination
if ( ! is.null(xsi.elim) ){
xsi.elim <- data.frame( xsi.elim, xsi.elim.index )
xsi.elim <- xsi.elim[ ! duplicated( xsi.elim[,2] ), ]
xsi.elim <- xsi.elim[ order( xsi.elim[,2] ), ]
}
#--- clean xsi.constr
xsi1 <- xsi.constr$xsi.constraints
xsi.constr$intercept_included <- FALSE
ind <- grep("(Intercept", rownames(xsi1), fixed=TRUE)
if ( length(ind) > 0 ){
xsi1 <- xsi1[ - ind, ]
xsi.constr$xsi.constraints <- xsi1
xsi.constr$intercept_included <- TRUE
}
xsi1 <- xsi.constr$xsi.table
ind <- grep("(Intercept", paste(xsi1$parameter), fixed=TRUE)
if ( length(ind) > 0 ){
xsi1 <- xsi1[ - ind, ]
xsi.constr$xsi.table <- xsi1
}
#cat(" +++ out .A.matrix" ) ; z1 <- Sys.time() ; print(z1-z0) ; z0 <- z1
res <- list( "A"=A, "X"=X.out, "otherFacets"=otherFacets, "xsi.constr"=xsi.constr,
"facet.design"=facet.design, "xsi.elim"=xsi.elim )
return(res)
}
# .A.matrix2 <- tam_A_matrix2
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.