R/tam.mml.mfr.R

In TAM: Test Analysis Modules

## File Name: tam.mml.mfr.R
## File Version: 9.955

tam.mml.mfr <- function( resp, Y=NULL, group=NULL,  irtmodel="1PL",
            formulaY=NULL, dataY=NULL, ndim=1, pid=NULL, xsi.fixed=NULL,
            xsi.setnull=NULL, xsi.inits=NULL, beta.fixed=NULL, beta.inits=NULL,
            variance.fixed=NULL, variance.inits=NULL, est.variance=TRUE,
            formulaA=~item+item:step, constraint="cases", A=NULL, B=NULL,
            B.fixed=NULL, Q=NULL, facets=NULL, est.slopegroups=NULL, E=NULL,
            pweights=NULL, verbose=TRUE, control=list(), delete.red.items=TRUE )
{
    CALL <- match.call()
    a0 <- Sys.time()
    s1 <- Sys.time()

    prior_list_xsi=NULL
    mstep_intercept_method <- "R"

    # display
    disp <- "....................................................\n"
    increment.factor <- progress <- nodes <- snodes <- ridge <- xsi.start0 <- QMC <- NULL
    maxiter <- conv <- convD <- min.variance <- max.increment <- Msteps <- convM <- NULL
    resp_orig <- resp
    B00 <- B
    B <- trim_increment <- NULL
    fac.oldxsi <- acceleration <- NULL

    #**** handle verbose argument
    args_CALL <- as.list( sys.call() )
    if ( ! tam_in_names_list( list=control, variable="progress" )     ){
        control$progress <- verbose
    }

    #--- attach control elements
    e1 <- environment()
    tam_fct <- "tam.mml.mfr"
    res <- tam_mml_control_list_define(control=control, envir=e1, tam_fct=tam_fct,
                prior_list_xsi=prior_list_xsi)
    con <- res$con
    con1a <- res$con1a

    #- check constraint
    constraint <- tam_mml_constraint_check(constraint=constraint)

    # userfct.variance is not allowed in tam.mml.mfr
    userfct.variance <- NULL

    fac.oldxsi <- max( 0, min( c( fac.oldxsi, .95 ) ) )
    if ( constraint=="items" ){ beta.fixed <- FALSE }

    pid0 <- pid <- unname(c(unlist(pid)))
    if (progress){
        cat(disp)
        cat("Processing Data     ", paste(Sys.time()), "\n") ; utils::flush.console()
    }
    if ( ! is.null(group) ){
      con1a$QMC <- QMC <- FALSE
      con1a$snodes <- snodes <- 0
    }

    resp <- as.matrix(resp)
    resp <- add.colnames.resp(resp)
    itemnames <- colnames(resp)

    nullY <- is.null(Y)

    if ( ! is.null(facets) ){
        facets <- as.data.frame(facets)
    }
# cat("read data" ) ; a1 <- Sys.time() ; print(a1-a0) ; a0 <- a1

    #--- compute maxKi
    res <- tam_mml_mfr_proc_compute_maxKi(resp=resp, facets=facets)
    maxKi <- res$maxKi

    #--- handle formula and facets
    resp00 <- resp
    res <- tam_mml_mfr_dataprep( formulaA=formulaA, xsi.setnull=xsi.setnull, B=B,
                Q=Q, resp=resp, pid=pid, facets=facets, beta.fixed=beta.fixed )
    formulaA <- res$formula_update
    xsi.setnull <- res$xsi.setnull
    beta.fixed <- res$beta.fixed
    facets <- res$facets
    PSF <- res$PSF
    pid <- res$pid
#cat(" mml mfr dataprep  " ) ; a1 <- Sys.time() ; print(a1-a0) ; a0 <- a1

    #--- create design matrices
    res <-  tam_mml_mfr_proc_create_design_matrices( pid=pid, maxKi=maxKi, resp=resp,
                formulaA=formulaA, facets=facets, constraint=constraint, ndim=ndim, Q=Q,
                A=A, B=B, progress=progress, xsi.fixed=xsi.fixed, resp00=resp00, B00=B00,
                beta.fixed=beta.fixed )
    pid <- res$pid
    diffKi <- res$diffKi
    var_ki <- res$var_ki
    xsi.fixed <- res$xsi.fixed
    xsi.elim <- res$xsi.elim
    beta.fixed <- res$beta.fixed
    A <- res$A
    cA <- res$cA
    B <- res$B
    Q <- res$Q
    X <- res$X
    X.red <- res$X.red
    gresp <- res$gresp
    gresp.noStep <- res$gresp.noStep
    xsi.constr <- res$xsi.constr
    design <- res$design
# cat(" --- design matrix ready" ) ; a1 <- Sys.time() ; print(a1-a0) ; a0 <- a1

    #--- processing in case of multiple person IDs in a dataset
    tp <- max(table(pid))
    if ( tp > 1){
        res <- tam_mml_mfr_proc_multiple_person_ids( pid=pid, tp=tp, gresp=gresp,
                    gresp.noStep=gresp.noStep, progress=progress, group=group, Y=Y,
                    pweights=pweights)
        pid <- res$pid
        gresp <- res$gresp
        gresp.noStep <- res$gresp.noStep
        group <- res$group
        Y <- res$Y
        pweights <- res$pweights
    }
# cat("process data in case of multiple persons" ) ; a1 <- Sys.time() ; print(a1-a0) ; a0 <- a1

    #--- set some xsi effects to zero
    res <- tam_mml_mfr_proc_xsi_setnull( xsi.setnull=xsi.setnull, A=A, xsi.fixed=xsi.fixed )
    xsi.fixed <- res$xsi.fixed
    xsi0 <- res$xsi0

    nitems <- nrow(X.red)
    nstud <- nrow(gresp)        # number of students
    if ( is.null(pweights) ){
        pweights <- rep(1,nstud) # weights of response pattern
    }

    if (progress){
        cat("    * Response Data:", nstud, "Persons and ",
            ncol(gresp.noStep), "Generalized Items (", paste(Sys.time()),")\n" )  ;
        utils::flush.console()
    }
    if ( is.null(pid) ){
        pid <- seq(1,nstud)
    }

    # normalize person weights to sum up to nstud
    pweights <- nstud * pweights / sum(pweights)
    # a matrix version of person weights
    pweightsM <- outer( pweights, rep(1,nitems) )

    # calculate ndim if only B or Q are supplied
    if ( ! is.null(B) ){ ndim <- dim(B)[3] }
    if ( ! is.null(Q) ){ ndim <- dim(Q)[2] }

    betaConv <- FALSE         #flag of regression coefficient convergence
    varConv <- FALSE          #flag of variance convergence
    nnodes <- length(nodes)^ndim
    if ( snodes > 0 ){ nnodes <- snodes }

    #--- print information about nodes
    res <- tam_mml_progress_proc_nodes( progress=progress, snodes=snodes, nnodes=nnodes,
                    skillspace="normal", QMC=QMC)

    #--- maximum no. of categories per item. Assuming dichotomous
    maxK <- max( resp, na.rm=TRUE ) + 1

    #--- number of parameters
    np <- dim(A)[[3]]

    #--- xsi parameter index
    res <- tam_mml_proc_est_xsi_index(A, xsi.inits, xsi.fixed)
    np <- res$np
    xsi <- res$xsi
    est.xsi.index0 <- est.xsi.index <- res$est.xsi.index

    #--- inits variance
    res <- tam_mml_inits_variance( variance.inits=variance.inits, ndim=ndim, variance.fixed=variance.fixed )
    variance <- res$variance

    #--- inits group
    res <- tam_mml_inits_groups( group=group )
    G <- res$G
    groups <- res$groups
    group <- res$group
    var.indices <- res$var.indices

    #--- inits beta
    res <- tam_mml_mfr_inits_beta( Y=Y, formulaY=formulaY, dataY=dataY, G=G, group=group, groups=groups,
                nstud=nstud, pweights=pweights, ridge=ridge, beta.fixed=beta.fixed, xsi.fixed=xsi.fixed,
                constraint=constraint, ndim=ndim, beta.inits=beta.inits, tp=tp, gresp=gresp,
                pid0=pid0 )
    Y <- res$Y
    nullY <- res$nullY
    formulaY <- res$formulaY
    nreg <- res$nreg
    W <- res$W
    YYinv <- res$YYinv
    beta.fixed <- res$beta.fixed
    beta <- res$beta

    #--- response indicators
    res <- tam_mml_mfr_proc_response_indicators(nitems, gresp, gresp.noStep)
    resp.ind.list <- res$resp.ind.list
    gresp.ind <- res$gresp.ind
    gresp.noStep.ind <- res$gresp.noStep.ind
    resp.ind <- res$resp.ind
    nomiss <- res$nomiss
    miss.items <- res$miss.items
    gresp0.noStep <- res$gresp0.noStep
    gresp <- res$gresp
    gresp.noStep <- res$gresp.noStep

    #-- delete items with only missing responses
    res <- tam_mml_mfr_proc_delete_missing_items( miss.items=miss.items,
                delete.red.items=delete.red.items, maxK=maxK,
                gresp=gresp, gresp.noStep=gresp.noStep, gresp.noStep.ind=gresp.noStep.ind,
                A=A, B=B, resp.ind.list=resp.ind.list, resp.ind=resp.ind, nitems=nitems,
                pweightsM=pweightsM, pweights=pweights, nstud=nstud, progress=progress )
    miss.itemsK  <- res$miss.itemsK
    miss.items <- res$miss.items
    delete.red.items <- res$delete.red.items
    A <- res$A
    B <- res$B
    gresp <- res$gresp
    gresp.noStep <- res$gresp.noStep
    gresp.noStep.ind <- res$gresp.noStep.ind
    resp.ind.list <- res$resp.ind.list
    resp.ind <- res$resp.ind
    nitems <- res$nitems
    pweightsM <- res$pweightsM

    #-- AXsi
    AXsi <- matrix(0,nrow=nitems,ncol=maxK )  #A times xsi

    #--- parameter indices xsi parameters
    res <- tam_mml_proc_xsi_parameter_index_A(A=A, np=np)
    indexIP <- res$indexIP
    indexIP.list <- res$indexIP.list
    indexIP.list2 <- res$indexIP.list2
    indexIP.no <- res$indexIP.no

    #--- sufficient statistics for item parameters
    cA <- t( matrix( aperm( A, c(2,1,3) ), nrow=dim(A)[3], byrow=TRUE ) )
    res <- tam_mml_sufficient_statistics( nitems=nitems, maxK=maxK, resp=gresp.noStep,
                    resp.ind=gresp.noStep.ind,     pweights=pweights, cA=cA, progress=progress )
    ItemScore <- res$ItemScore
    cResp <- res$cResp
    col.index <- res$col.index

    #--- inits xsi
    res <- tam_mml_mfr_inits_xsi( gresp.noStep.ind=gresp.noStep.ind, col.index=col.index, cA=cA,
                pweights=pweights, xsi=xsi, xsi.start0=xsi.start0, resp=resp, A=A,
                xsi.inits=xsi.inits, xsi.fixed=xsi.fixed, ItemScore=ItemScore, est.xsi.index=est.xsi.index )
    xsi <- res$xsi
    ItemMax <- res$ItemMax

    #--- prior distribution xsi
    prior_list_xsi <- tam_mml_proc_prior_list_xsi( prior_list_xsi=prior_list_xsi, xsi=xsi )

    xsi.min.deviance <- xsi
    beta.min.deviance <- beta
    variance.min.deviance <- variance

    #--- create grid of nodes for numeric or stochastic integration
    res <- tam_mml_create_nodes( snodes=snodes, nodes=nodes, ndim=ndim, QMC=QMC )
    theta <- res$theta
    theta2 <- res$theta2
    thetawidth <- res$thetawidth
    theta0.samp <- res$theta0.samp
    thetasamp.density <- res$thetasamp.density

    deviance <- 0
    deviance.history <- tam_deviance_history_init(maxiter=maxiter)

    iter <- 0
    a02 <- a1 <- 999    # item parameter change
    a4 <- 0

    hwt.min <- 0
    rprobs.min <- 0
    AXsi.min <- 0
    B.min <- 0
    deviance.min <- 1E100
    itemwt.min <- 0

    #--- create unidim_simplify
    res <- tam_mml_proc_unidim_simplify( Y=Y, A=A, G=G, beta.fixed=beta.fixed )
    unidim_simplify <- res$unidim_simplify
    YSD <- res$YSD
    Avector <- res$Avector

    #--- acceleration
    res <- tam_acceleration_inits(acceleration=acceleration, G=G, xsi=xsi,
                variance=variance)
    xsi_acceleration <- res$xsi_acceleration
    variance_acceleration <- res$variance_acceleration

    #--- warning multiple group estimation
    res <- tam_mml_warning_message_multiple_group_models( ndim=ndim, G=G)

    #--- compute some arguments for EM algorithm
    maxcat <- tam_rcpp_mml_maxcat(A=as.vector(A), dimA=dim(A) )

    ##**SE
    se.xsi <- 0*xsi
    se.B <- 0*B
    se.xsi.min <- se.xsi
    se.B.min <- se.B
    devch <- 0

    # display
    disp <- "....................................................\n"
    # define progress bar for M step
# cat("rest  " ) ; a1 <- Sys.time() ; print(a1-a0) ; a0 <- a1

    ##############################################################
    ##############################################################
    ##############################################################
    #Start EM loop here
    while ( ( (!betaConv | !varConv)  | ((a1 > conv) | (a4 > conv) | (a02 > convD)) )  &
              (iter < maxiter) ) {

        # a0 <- Sys.time()
        iter <- iter + 1
        #--- progress
        res <- tam_mml_progress_em0(progress=progress, iter=iter, disp=disp)
        # calculate nodes for Monte Carlo integration
        if ( snodes > 0){
            res <- tam_mml_update_stochastic_nodes( theta0.samp=theta0.samp, variance=variance,
                        snodes=snodes, beta=beta, theta=theta )
            theta <- res$theta
            theta2 <- res$theta2
            thetasamp.density <- res$thetasamp.density
        }
        olddeviance <- deviance
        #--- calculation of probabilities
        res <- tam_mml_calc_prob( iIndex=1:nitems, A=A, AXsi=AXsi, B=B, xsi=xsi, theta=theta, nnodes=nnodes,
                    maxK=maxK, recalc=TRUE, maxcat=maxcat, use_rcpp=TRUE )
      # cat("calc prob") ; a1 <- Sys.time(); print(a1-a0) ; a0 <- a1
        rprobs <- res$rprobs
        AXsi <- res$AXsi

        #--- calculate student's prior distribution
        gwt <- tam_stud_prior( theta=theta, Y=Y, beta=beta, variance=variance, nstud=nstud,
                    nnodes=nnodes, ndim=ndim, YSD=YSD, unidim_simplify=unidim_simplify,
                    snodes=snodes )

        #--- calculate student's likelihood
        res.hwt <- tam_calc_posterior( rprobs=rprobs, gwt=gwt, resp=gresp.noStep, nitems=nitems,
                        resp.ind.list=resp.ind.list, normalization=TRUE,
                        thetasamp.density=thetasamp.density, snodes=snodes, resp.ind=resp.ind,
                        avoid.zerosum=TRUE )
        hwt <- res.hwt$hwt

        #--- M step: estimation of beta and variance
        resr <- tam_mml_mstep_regression( resp=gresp.noStep, hwt=hwt,
                    resp.ind=gresp.noStep.ind, pweights=pweights, pweightsM=pweightsM,
                    Y=Y, theta=theta, theta2=theta2, YYinv=YYinv, ndim=ndim, nstud=nstud,
                    beta.fixed=beta.fixed, variance=variance, Variance.fixed=variance.fixed,
                    group=group, G=G, snodes=snodes, nomiss=nomiss, iter=iter,
                    min.variance=min.variance, userfct.variance=userfct.variance,
                    variance_acceleration=variance_acceleration, est.variance=est.variance,
                    beta=beta )
        beta <- resr$beta
        variance <- resr$variance
        itemwt <- resr$itemwt
        variance_acceleration <- resr$variance_acceleration
        variance_change <- resr$variance_change
        beta_change <- resr$beta_change

        if ( beta_change < conv){ betaConv <- TRUE }
        if ( variance_change < conv){ varConv <- TRUE }

        #--- M-step item intercepts
        if (mstep_intercept_method=="optim"){
            res <- tam_calc_counts( resp=gresp.noStep, theta=theta, resp.ind=gresp.noStep.ind, group=group,
                        maxK=maxK, pweights=pweights, hwt=hwt )
            n.ik <- res$n.ik
        }
        res <- tam_mml_mstep_intercept( A=A, xsi=xsi, AXsi=AXsi, B=B, theta=theta,
                    nnodes=nnodes, maxK=maxK, Msteps=Msteps, rprobs=rprobs, np=np,
                    est.xsi.index0=est.xsi.index0, itemwt=itemwt, indexIP.no=indexIP.no,
                    indexIP.list2=indexIP.list2, Avector=Avector, max.increment=max.increment,
                    xsi.fixed=xsi.fixed, fac.oldxsi=fac.oldxsi, ItemScore=ItemScore,
                    convM=convM, progress=progress, nitems=nitems, iter=iter,
                    increment.factor=increment.factor, xsi_acceleration=xsi_acceleration,
                    trim_increment=trim_increment, prior_list_xsi=prior_list_xsi,
                    mstep_intercept_method=mstep_intercept_method, n.ik=n.ik, maxcat=maxcat )
        xsi <- res$xsi
        se.xsi <- res$se.xsi
        max.increment <- res$max.increment
        xsi_acceleration <- res$xsi_acceleration
        xsi_change <- res$xsi_change
        logprior_xsi <- res$logprior_xsi

        #--- compute deviance
        res <- tam_mml_compute_deviance( loglike_num=res.hwt$rfx, loglike_sto=res.hwt$rfx,
                    snodes=snodes, thetawidth=thetawidth, pweights=pweights, deviance=deviance,
                    deviance.history=deviance.history, iter=iter, logprior_xsi=logprior_xsi )
        deviance <- res$deviance
        deviance.history <- res$deviance.history
        a01 <- rel_deviance_change <- res$rel_deviance_change
        a02 <- deviance_change <- res$deviance_change
        if (con$dev_crit=="relative" ){ a02 <- a01 }
        penalty_xsi <- res$penalty_xsi
        deviance_change_signed <- res$deviance_change_signed

        if( deviance < deviance.min ){
            xsi.min.deviance <- xsi
            beta.min.deviance <- beta
            variance.min.deviance <- variance
            hwt.min <- hwt
            AXsi.min <- AXsi
            B.min <- B
            deviance.min <- deviance
            itemwt.min <- itemwt
            se.xsi.min <- se.xsi
            se.B.min <- se.B
        }

        a1 <- xsi_change
        a2 <- beta_change
        a3 <- variance_change
        devch <- - ( deviance - olddeviance )

        #--- print progress
        res <- tam_mml_progress_em( progress=progress, deviance=deviance, deviance_change=deviance_change,
                    iter=iter, rel_deviance_change=rel_deviance_change, xsi_change=xsi_change,
                    beta_change=beta_change, variance_change=variance_change, B_change=0,
                    devch=devch, penalty_xsi=penalty_xsi )

    } # end of EM loop
    #############################################################
    #############################################################
    xsi.min.deviance -> xsi
    beta.min.deviance -> beta
    variance.min.deviance -> variance
    hwt.min -> hwt
    AXsi.min -> AXsi
    B.min -> B
    deviance.min -> deviance
    itemwt.min -> itemwt
    se.xsi.min -> se.xsi
    se.B.min -> se.B
    #******
    #***
    resp <- gresp0.noStep
    resp.ind <- gresp.noStep.ind

    #*** include NAs in AXsi
    AXsi <- tam_mml_include_NA_AXsi(AXsi=AXsi, maxcat=maxcat, A=A, xsi=xsi)

    #****
    # look for non-estimable xsi parameters
#    xsi[ xsi==99 ] <- NA

    #******
    # generate input for fixed parameters
    xsi.fixed.estimated <- tam_generate_xsi_fixed_estimated( xsi=xsi, A=A )
    B.fixed.estimated <- tam_generate_B_fixed_estimated(B=B)

    #**** standard errors AXsi
    se.AXsi <- tam_mml_se_AXsi( AXsi=AXsi, A=A, se.xsi=se.xsi, maxK=maxK )

    ##*** information criteria
    ic <- tam_mml_ic( nstud=nstud, deviance=deviance, xsi=xsi, xsi.fixed=xsi.fixed,
                beta=beta, beta.fixed=beta.fixed, ndim=ndim,
                variance.fixed=variance.fixed, G=G, irtmodel=irtmodel, B_orig=NULL,
                B.fixed=B.fixed, E=E, est.variance=TRUE, resp=resp,
                est.slopegroups=NULL, variance.Npars=NULL, group=group,
                penalty_xsi=penalty_xsi, pweights=pweights, resp.ind=resp.ind )


    #*** calculate counts
    res <- tam_calc_counts( resp=gresp.noStep, theta=theta, resp.ind=gresp.noStep.ind,
                group=group, maxK=maxK, pweights=pweights, hwt=hwt )
    n.ik <- res$n.ik
    pi.k <- res$pi.k

    #--- collect item parameters
    item1 <- tam_itempartable( resp=gresp.noStep, maxK=maxK, AXsi=AXsi, B=B,
                    ndim=ndim, resp.ind=gresp.noStep.ind,
                    rprobs=rprobs, n.ik=n.ik, pi.k=pi.k, order=TRUE, pweights=pweights )

    #--- collect all person statistics
    res <- tam_mml_person_posterior( pid=pid, nstud=nstud, pweights=pweights,
                resp=gresp.noStep, resp.ind=gresp.noStep.ind, snodes=snodes,
                hwtE=hwt, hwt=hwt, ndim=ndim, theta=theta )
    person <- res$person
    EAP.rel <- res$EAP.rel

    #******
    s2 <- Sys.time()
    item <- data.frame( "xsi.index"=1:np, "xsi.label"=dimnames(A)[[3]], "est"=xsi )
    if (progress){
        cat(disp)
        cat("Item Parameters\n")
        item2 <- item
        item2[,"est"] <- round( item2[,"est"], 4 )
        print(item2)
        cat("...................................\n")
        cat("Regression Coefficients\n")
        print( beta, 4  )
        cat("\nVariance:\n" ) #, round( varianceM, 4 ))
        if (G==1 ){
            varianceM <- matrix( variance, nrow=ndim, ncol=ndim )
            print( varianceM, 4 )
        } else {
            print( variance[ var.indices], 4 )
        }
        if ( ndim > 1){
                cat("\nCorrelation Matrix:\n" ) #, round( varianceM, 4 ))
            print( stats::cov2cor(varianceM), 4 )
        }
        cat("\n\nEAP Reliability:\n")
        print( round (EAP.rel,3) )
        cat("\n-----------------------------")
        devmin <- which.min( deviance.history[,2] )
        if ( devmin < iter ){
            cat(paste("\n\nMinimal deviance at iteration ", devmin,
                  " with deviance ", round(deviance.history[ devmin, 2 ],3), sep=""), "\n")
            cat("The corresponding estimates are\n")
            cat("  xsi.min.deviance\n  beta.min.deviance \n  variance.min.deviance\n\n")
        }
        cat( "\nStart: ", paste(s1))
        cat( "\nEnd: ", paste(s2),"\n")
        print(s2-s1)
        cat( "\n" )
    }

    #--- collect xsi parameters
    res <- tam_mml_mfr_collect_xsi_parameters( xsi.constr=xsi.constr, resp=resp, A=A, xsi=xsi,
                se.xsi=se.xsi, delete.red.items=delete.red.items, itemnames=itemnames,
                miss.items=miss.items )
    resp <- res$resp
    xsi <- res$xsi
    xsi.facets <- res$xsi.facets

    #--- recompute posterior
    res.hwt <- tam_calc_posterior( rprobs=rprobs, gwt=1+0*gwt, resp=resp, nitems=nitems,
                    resp.ind.list=resp.ind.list, normalization=FALSE, thetasamp.density=thetasamp.density,
                    snodes=snodes, resp.ind=resp.ind )
    res.like <- res.hwt$hwt

    #***** standardized coefficients
    latreg_stand <- tam_latent_regression_standardized_solution(variance=variance, beta=beta, Y=Y)

    #--- OUTPUT LIST
    deviance.history <- deviance.history[ 1:iter, ]
    res <- list( "xsi"=xsi, "xsi.facets"=xsi.facets,
                 "beta"=beta, "variance"=variance,
                 "item"=item1,
                 "person"=person, pid=pid, "EAP.rel"=EAP.rel,
                 "post"=hwt,  "rprobs"=rprobs, "itemweight"=itemwt,
                 "theta"=theta,
                 "n.ik"=n.ik, "pi.k"=pi.k,
                 "Y"=Y, "resp"=resp,
                 "resp.ind"=resp.ind, "group"=group,
                 "G"=if ( is.null(group)){1} else { length(unique( group ) )},
                 "groups"=if ( is.null(group)){1} else { groups },
                 "formulaY"=formulaY, "dataY"=dataY,
                 "pweights"=pweights,
                 "time"=c(s1,s2), "A"=A, "B"=B,
                 "se.B"=se.B,
                 "nitems"=nitems, "maxK"=maxK, "AXsi"=AXsi,
                 "AXsi_"=- AXsi,
                 "se.AXsi"=se.AXsi,
                 "nstud"=nstud, "resp.ind.list"=resp.ind.list,
                 "hwt"=hwt, "like"=res.like, "ndim"=ndim,
                 "xsi.fixed"=xsi.fixed,
                 "xsi.fixed.estimated"=xsi.fixed.estimated,
                 "B.fixed.estimated"=B.fixed.estimated,
                 "beta.fixed"=beta.fixed, "Q"=Q,
                 "formulaA"=formulaA, "facets"=facets,
                 "xsi.constr"=xsi.constr,
                 "variance.fixed"=variance.fixed,
                 "nnodes"=nnodes, "deviance"=ic$deviance,
                 "ic"=ic, thetasamp.density=thetasamp.density,
                 "deviance.history"=deviance.history,
                 "control"=con1a, "irtmodel"=irtmodel,
                 "iter"=iter, "resp_orig"=resp_orig,
                 "printxsi"=TRUE, "YSD"=YSD, "PSF"=PSF,
                 CALL=CALL, latreg_stand=latreg_stand,
                  prior_list_xsi=prior_list_xsi, penalty_xsi=penalty_xsi
    )
    class(res) <- "tam.mml"
    return(res)
 }