Nothing
### Roxygen-izable Documentation ----
#' Create model matrices for Thurstonian modeling
#'
#' This function creates and returns model matrices of Thurstonian model equations.
#'
#'
### Inputs ----
#'
#' @param blocks A matrix defining the blocks of the model. The number of rows
#' must be the number of blocks, each row represents a block and contains the
#' item numbers. The number of columns present the number of items per block.
#' @param itf A vector defining the items-to-factor relation. For example
#' `c(1,1,1,2,2,2)` defines six items, the first three correspond to factor 1,
#' the second three correspond to factor 2.
#' @param model A descriptor for the model. Can be one of `'lmean'`,
#' `'uc'`, `'irt'` or `'simple2'`, `'simple3'` or `'simple5'`. The Number behind
#' the `'simple'` statement defines the Thurstone case.
#'
### Outputs ----
#' @return Returns a list of elements containing model matrix information.
#'
### Examples ----
#' @examples
#'
#' # set seed and define blocks
#' set.seed(1)
#' blocks <- matrix(sample(1:15,15), ncol = 3)
#'
#' # define the item-to-factor relation
#' itf <- rep(1:3,5)
#'
#' mod.matrices(blocks,itf,'irt')
#'
#' @export
### Function definition ----
mod.matrices <- function(blocks,itf,model){
if(model=='simple2'){
model <- 'simple'
case <- 2
}
if(model=='simple3'){
model <- 'simple'
case <- 3
}
if(model=='simple5'){
model <- 'simple'
case <- 5
}
if(!is.matrix(blocks)&!is.vector(blocks)){
stop('Blocks must be given as a matrix. Only for full designs a vector is possible.')
}
if(is.vector(blocks)){
blocks <- blocks <- matrix(blocks,nrow=1)
}
design <- ifelse(dim(blocks)[1]==1,'full','block')
# Find number of items per block
nrank <- dim(blocks)[2]
# Find number of factors, blocks and items
nblock <- dim(blocks)[1]
nitem <- length(unique(c(blocks)))
nfactor <- max(itf)
# Pair combinations
pair_combn <- t(utils::combn(1:nitem,2))
# Define items to factors for lambda (loadings)
lambda_raw <- matrix(0,nrow=nitem,ncol=nfactor)
for(i in 1:nfactor){
lambda_raw[,i] <- ifelse(itf==i,i,0)
}
# Create names of pairs
pair_names <- paste0('i',pair_combn[,1],'i',pair_combn[,2])
# For each block: is it uni-dimensional?
blocks_lambda <- matrix(itf[blocks],ncol=nrank)
unidim <- apply(blocks_lambda,1,function(x) min(x)==max(x))
# Create names of pairs for blocks
pair_names_bs <- i.name(blocksort(blocks))
pair_names_b <- i.name(blocks)
tmp <- which(pair.combn(blocks)[,1]>pair.combn(blocks)[,2])
if(length(tmp)!=0){
tmp1 <- pair_names_b[tmp]
tmp2 <- sub('^i.+i','i',tmp1)
tmp3 <- tmp1
for(j in 1:length(tmp)){
tmp3 <- paste0(tmp2[j],sub(paste0(tmp2[j],'$'),'',tmp1[j]))
pair_names_bs[which(pair_names_bs%in%tmp3)] <- pair_names_b[tmp[j]]
}
}
pair_combn_b <- pair.combn(blocks)
# MODEL
if(model!='simple'){
traits <- paste0('Trait',1:nfactor)
}
# Create Matrix A
A <- designA(blocks)
Arref <- rref(A)
Arank <- dim(A)[1]-length(which(apply(abs(Arref),1,sum)==0))
# Find global blocks
gblocks <- metablock(blocks)
unidim_gblocks <- sapply(gblocks,function(x) min(itf[x])==max(itf[x]))
# Latent utility model and the unconstrained threshold model
if(model%in%c('lmean','uc','simple')){
# Part 1 (primary factors)
t <- paste0('t',1:nitem)
A_sign <- ifelse(A==0,'',ifelse(A==1,'+','-'))
At <- matrix(NA, nrow=nrow(A_sign), ncol=nitem)
At <- t(apply(A_sign,1,function(x)sub('\\+','',sub('^_0.*','',sub('^\\+0.*','',sub('^0.*','',paste0(x,t)))))))
At <- sub('^-','',matrix(paste0(At,' BY ',pair_names_b,'@',ifelse(A_sign=='+','',A_sign),'1;\n'), ncol=nitem))
At <- ifelse(A_sign=='','',At)
At_text <- paste0(apply(At,2,function(x)paste0(x,collapse = '')), collapse = '')
# Part 2 (secondary factors)
if(model!='simple'){
if(nfactor==1&Arank==(nitem-1)){
input <- paste0(input, 'Trait1 BY\nt1@1\n', paste0(paste0('t',2:nitem,'*1'),collapse='\n'), ';\n\n')
}else{
A_fac <- rep('',nfactor)
itf_gblocks <- lapply(gblocks,function(x) itf[x])
items_fac <- list()
items_fac_uni <- list()
for(i in 1:nfactor){
item_fac <- c()
item_fac_uni <- c()
for(j in 1:length(gblocks)){
if(unidim_gblocks[j]){
item_fac_uni <- c(item_fac_uni,gblocks[[j]][which(itf_gblocks[[j]]==i)[1]])
item_fac <- c(item_fac,gblocks[[j]][which(itf_gblocks[[j]]==i)[-1]])
}else{
item_fac <- c(item_fac,gblocks[[j]][which(itf_gblocks[[j]]==i)])
}
}
items_fac[[i]] <- sort(item_fac)
if(all(is.na(item_fac_uni))){
items_fac_uni[[i]] <- NA
}else{
items_fac_uni[[i]] <- sort(item_fac_uni)
}
A_fac[[i]] <- paste0(traits[i], ' BY\n')
A_fac[[i]] <- paste0(paste0(traits[i], ' BY\n'),paste0('t',unname(items_fac[[i]]),'*1', '\n', collapse = ''))
#for(j in 1:length(gblocks)){
#if(unidim_gblocks[j]){
#if(length(items_fac_uni[[i]])>0){
if(all(!is.na(items_fac_uni[[i]]))){
A_fac[[i]] <- paste0(A_fac[[i]],paste0('t',unname(items_fac_uni[[i]]),'@1', '\n', collapse = ''))
}
#}
#}
#}
}
}
}
}
# IRT model
if(model%in%c('irt')){
# Create matrix lambda
lambda <- ifelse(lambda_raw==0,0,paste0('L',row(lambda_raw)))
A_sign <- ifelse(A==0,0,ifelse(A==1,'+','_'))
A_lambda <- A%*%ifelse(lambda_raw==0,0,1)
Lambda <- matrix(NA, nrow=nrow(A_sign), ncol=ncol(lambda))
for(i in 1:nfactor){
Lambda[,i] <- apply(A_sign,1,function(x) paste0(sub('^_0.*','',sub('^\\+0.*','',sub('^0.*','',paste0(x,lambda[,i])))), collapse =''))
}
tmp <- grep('^_L.+\\+',Lambda)
tmp1 <- Lambda[tmp]
Lambda[tmp] <- paste0(sub('^_L.+\\+','',tmp1),sub('\\+L.+','',tmp1))
Lambda <- sub('\\+','',Lambda)
Lambda <- sub('^_L','L_',Lambda)
A_mplus <- NULL
for(i in 1:nfactor){
A_mplus <- paste0(A_mplus,'\n\n',traits[i],' BY \n',
paste0(ifelse(Lambda[,i]=='','',
paste0(pair_names_b,'*',A_lambda[,i],' (',Lambda[,i],')\n')),collapse = ''),';')
}
if(nrank==2){
Psi <- ifelse(row(diag(nitem))==col(diag(nitem)),paste0('1/2'),0)
Psi <- ifelse(Psi==0,'',Psi) # neu
A_sign <- ifelse(A==0,0,ifelse(A==1,'+','-')) # alt
APsi <- matrix(NA, nrow=length(pair_names_b), ncol=nitem)
for(i in 1:nitem){
for(j in 1:length(pair_names_b)){
tmp <- paste0(A_sign[j,],Psi[,i])
tmp <- sub('^-$','',sub('^\\+$','',tmp))
tmp <- sub('^-0.*','',sub('^\\+0.*','',sub('^0.*','',tmp)))
APsi[j,i] <- paste0(tmp, collapse ='')
}
}
APsi <- sub('^\\+','',APsi)
APsiA <- matrix(NA, nrow=length(pair_names_b), ncol=length(pair_names_b))
for(i in 1:length(pair_names_b)){
for(j in 1:length(pair_names_b)){
tmp <- paste0(A_sign[i,],APsi[j,])
#tmp <- sub('\\+','',sub('^0.*','',tmp))
tmp <- sub('^0.*','',tmp)
APsiA[j,i] <- paste0(tmp, collapse ='')
}
}
APsiA <- ifelse(APsiA=='+-','',ifelse(APsiA=='-+','',APsiA))
APsiA <- gsub('-\\+','',APsiA)
APsiA <- gsub('\\+-','-',APsiA)
APsiA <- gsub('--','+',APsiA)
APsiA <- sub('-$','',APsiA)#
APsiA <- sub('\\+$','',APsiA)#
APsiA <- sub('^\\+','',APsiA)
APsiA <- sub('1/2+1/2','1',APsiA,fixed = TRUE)
if(length(gblocks)!=nblock){
# Define new parameters
unidim_item <- Lambda[grep('^L.*_L',Lambda)]
tmp <- unique(c(sub('_.*','',unidim_item),sub('^.*_','',unidim_item)))
if(length(which(tmp%in%Lambda))!=0){
tmp <- tmp[-which(tmp%in%Lambda)]
}
AddL <- c(paste0('L',1:nitem)[which(!paste0('L',1:nitem)%in%Lambda)])
tmp <- c(AddL,tmp)
Add_param <- unique(c(sort(tmp[which(nchar(tmp)==2)]),sort(tmp[which(nchar(tmp)==3)]),sort(tmp[which(nchar(tmp)==4)])))
Add_unidim <- NULL
for(i in 1:nrow(blocks)){
tmp <- blocks[i,]%in%sub('L','',Add_param)
if(all(tmp)){
Add_unidim <- c(Add_unidim ,paste0('L',blocks[i,]))
}
}
# Constraints
conL <- rep('',nitem)
for(i in 1:nitem){
if(paste0('L',i)%in%Lambda|paste0('L',i)%in%Add_param){
if(paste0('L_',i)%in%Lambda){
conL[i] <- paste0('L_',i,' = ','- L',i)
}
}else if(i%in%sub('L','',Add_param[-which(Add_param%in%Add_unidim)])){
conL[i] <- paste0('L_',i,' = ','- L',i)
}
}
Add_param <- paste0(paste0(Add_param, collapse = '\n'),';')
unidim_item <- ifelse(length(unidim_item)==0,'',paste0(unidim_item,' = ', sub('_','-',unidim_item),';\n',collapse = ''))
conL <- paste0(ifelse(conL=='','',paste0(conL,';\n')), collapse = '')
}else{
Add_param <- NULL
}
}
if(nrank>2){
# Create Matrix Psi for irt models
Psi <- ifelse(row(diag(nitem))==col(diag(nitem)),paste0('e',1:nitem),0)
Psi <- ifelse(Psi==0,'',Psi) # neu
A_sign <- ifelse(A==0,0,ifelse(A==1,'+','-')) # alt
APsi <- matrix(NA, nrow=length(pair_names_b), ncol=nitem)
for(i in 1:nitem){
for(j in 1:length(pair_names_b)){
tmp <- paste0(A_sign[j,],Psi[,i])
tmp <- sub('^-$','',sub('^\\+$','',tmp))
tmp <- sub('^-0.*','',sub('^\\+0.*','',sub('^0.*','',tmp)))
APsi[j,i] <- paste0(tmp, collapse ='')
}
}
APsi <- sub('^\\+','',APsi)
APsiA <- matrix(NA, nrow=length(pair_names_b), ncol=length(pair_names_b))
for(i in 1:length(pair_names_b)){
for(j in 1:length(pair_names_b)){
tmp <- paste0(A_sign[i,],APsi[j,])
tmp <- sub('^0.*','',tmp)
APsiA[j,i] <- paste0(tmp, collapse ='')
}
}
APsiA <- ifelse(APsiA=='+-','',ifelse(APsiA=='-+','',APsiA))
APsiA <- gsub('-\\+','',APsiA)
APsiA <- gsub('\\+-','-',APsiA)
APsiA <- gsub('--','+',APsiA)
APsiA <- sub('-$','',APsiA)#
APsiA <- sub('\\+$','',APsiA)#
APsiA <- sub('^\\+','',APsiA)
APsiA_under <- sub('-e','e_',APsiA)
APsiA_sign <- ifelse(APsiA=='','',ifelse(substr(APsiA,1,1)=='-','-',''))
# Define new parameters
unidim_item <- Lambda[grep('^L.*_L',Lambda)]
tmp <- unique(c(sub('_.*','',unidim_item),sub('^.*_','',unidim_item)))
if(length(which(tmp%in%Lambda))!=0){
tmp <- tmp[-which(tmp%in%Lambda)]
}
AddL <- c(paste0('L',1:nitem)[which(!paste0('L',1:nitem)%in%Lambda)])
tmp <- c(AddL,tmp)
Add_param <- unique(c(sort(tmp[which(nchar(tmp)==2)]),sort(tmp[which(nchar(tmp)==3)]),sort(tmp[which(nchar(tmp)==4)])))
Add_unidim <- NULL
for(i in 1:nrow(blocks)){
tmp <- blocks[i,]%in%sub('L','',Add_param)
if(all(tmp)){
Add_unidim <- c(Add_unidim ,paste0('L',blocks[i,]))
}
}
# Constraints
conL <- rep('',nitem)
for(i in 1:nitem){
if(paste0('L',i)%in%Lambda|paste0('L',i)%in%Add_param){
if(paste0('L_',i)%in%Lambda){
conL[i] <- paste0('L_',i,' = ','- L',i)
}
}else if(i%in%sub('L','',Add_param[-which(Add_param%in%Add_unidim)])){
conL[i] <- paste0('L_',i,' = ','- L',i)
}
}
Add_param <- paste0(paste0(Add_param, collapse = '\n'),';')
unidim_item <- ifelse(length(unidim_item)==0,'',paste0(unidim_item,' = ', sub('_','-',unidim_item),';\n',collapse = ''))
conL <- paste0(ifelse(conL=='','',paste0(conL,';\n')), collapse = '')
# Pair's uniqueness is equal to sum of two utility uniquenesses
if(design=='full'){
}else if(design=='block'){
uniquenesses <- c()
for(i in 1:nrow(blocks)){
items <- blocks[i,]
whitems <- which(pair_combn_b[,1]%in%items&pair_combn_b[,2]%in%items)
tmp <- pair_combn_b[whitems,1][which(pair_combn_b[whitems,1]%in%pair_combn_b[whitems,2])]
tmp <- paste0('e',pair_combn_b[whitems,1],'e',pair_combn_b[whitems,2],' = ' ,ifelse(pair_combn_b[whitems,1]%in%tmp,'-e_','e'),pair_combn_b[whitems,1],
ifelse(pair_combn_b[whitems,2]%in%tmp,' - e_',' + e'),pair_combn_b[whitems,2],';')
tmp0 <- paste0('e',pair_combn_b[whitems,1],'e',pair_combn_b[whitems,2])
for(i in length(tmp0):1){
if(any(grepl(tmp0[i],uniquenesses))){
tmp <- tmp[-i]
}
}
uniquenesses <- c(uniquenesses,tmp)
}
uniquenesses <- unique(uniquenesses[order(nchar(uniquenesses),uniquenesses)])
}
}
}
if(model=='irt'){
return(list(lambda=Lambda,
resvar=paste0(pair_names_b,'*2 ','(e',pair_combn_b[,1],'e',pair_combn_b[,2],');', collapse = '\n'),
Psi=Psi,
APsi=APsi,
APsiA=APsiA,
npair=length(pair_names_b),
naddparam=length(unlist(strsplit(gsub(';','',Add_param),'\n'))),
addparam=strsplit(gsub(';','',Add_param),'\n'),
pair_names_b=sub('i','I',pair_names_b)))
}else if(model=='lmean'|model=='uc'){
lambda_read <- strsplit(A_fac,'\n')
#if(connected|xtra){
# for(i in 1:nfactor){
# lambda_read[[i]] <- lambda_read[[i]][1:(length(which(itf==i))+1)]
# }
#}
return(list(lambda=lambda_read,
npair=length(pair_names_b),
pair_names_b=sub('i','I',pair_names_b)))
}else{
#lambda_read <- strsplit(A_fac,'\n')
#if(connected|xtra){
# for(i in 1:nfactor){
# lambda_read[[i]] <- lambda_read[[i]][1:(length(which(itf==i))+1)]
# }
#}
return(list(npair=length(pair_names_b),
pair_names_b=sub('i','I',pair_names_b)))
}
}
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.