version 0.8.6.4

DESCRIPTION

@@ -2,8 +2,8 @@ Package: mclogit
 Type: Package
 Title: Multinomial Logit Models, with or without Random Effects or
         Overdispersion
-Version: 0.8.5.1
-Date: 2020-06-27
+Version: 0.8.6.4
+Date: 2020-11-03
 Author: Martin Elff
 Maintainer: Martin Elff <mclogit@elff.eu>
 Description: Provides estimators for multinomial logit models in their
@@ -18,10 +18,10 @@ Imports: memisc, methods
 Suggests: MASS, nnet
 LazyLoad: Yes
 URL:
-        http://www.elff.eu/software/mclogit/,http://github.com/melff/mclogit/
-BugReports: http://github.com/melff/mclogit/issues
+        http://www.elff.eu/software/mclogit/,https://github.com/melff/mclogit/
+BugReports: https://github.com/melff/mclogit/issues
 RoxygenNote: 7.1.0
 NeedsCompilation: no
-Packaged: 2020-06-28 19:23:01 UTC; elff
+Packaged: 2020-11-04 11:23:11 UTC; elff
 Repository: CRAN
-Date/Publication: 2020-07-01 16:40:07 UTC
+Date/Publication: 2020-11-04 12:30:02 UTC

MD5

@@ -1,28 +1,29 @@
-17afe8ac573c522861b4839d631665b1 *DESCRIPTION
-4571d2fcfee83a95a2da51a09be7689e *NAMESPACE
+79c67f9b50f290802cee35b37ca7efde *DESCRIPTION
+253457c7632c75e465a2a23d604b0a49 *NAMESPACE
 5f46b25a049e81e86774e78cd294ff89 *R/AIC-mclogit.R
 a094e3ab8cf76cb83c81915e3199f5b3 *R/anova-mclogit.R
-6b3ff7bb0a6aa03542732e50a57a3936 *R/blockMatrices.R
+a807d113ca4268914fc1db4e4ef06d01 *R/blockMatrices.R
 825e92e98bee4d5cba8ebbb5bad2f448 *R/getSummary-mblogit.R
 172c1ad095181d15a3ac6267422d4e49 *R/getSummary-mclogit.R
-2326ed8abc8fe61028d6e72f58302dfe *R/mblogit.R
-2153774153109c8588f8743d3ca54bfe *R/mclogit-dispersion.R
-97c2d4ce949e3e1fe6ea6f79318e7770 *R/mclogit-fit.R
-73665a63241f8c02229a95bdad41e5cc *R/mclogit.R
-17ad276f4954104fe102a31abb666dce *R/mmclogit-fitPQLMQL.R
+fba3085b148d1da2e4fbcfdc53c4b191 *R/mblogit.R
+f814c3a2ac1bd0a7a0d953085bb68f4d *R/mclogit-dispersion.R
+b3caa0afcf4772d8b8e079626be1d497 *R/mclogit-fit.R
+8eea731951904e11a569c2cf05991ebd *R/mclogit.R
+cd9fc5997187b9a0919c7c9c14cbabfb *R/mmclogit-fitPQLMQL.R
 202e8dd6dd5ee19a281df6277484e46c *R/zzz.R
-fe1625a190f8a0e4ebec50342fcd3510 *data/Transport.rda
-8f626af2a1cc35f0632fec28e604d404 *data/electors.rda
+57ffbe9a81cbc936c04fa7e4d42c8e72 *data/Transport.rda
+f01c2450d12eccf5241c10a2040cd0b8 *data/electors.rda
 de4c63b49c62998527d72c185ceab0d3 *demo/00Index
 6ba2e10f92f3743f465cda6cc059f5a1 *demo/mclogit.test.R
-abeab02da03df6667eab7f296fcb4051 *inst/ChangeLog
-f1e22d308d70ae4b5403f5148132f255 *inst/NEWS.Rd
+e96435646fe124da3583fe9a2f005e2c *inst/ChangeLog
+2606b8800e9bf4d36effa9f523532c50 *inst/NEWS.Rd
 c05ff442f7a2312ff8f3d6b5aa96ddeb *man/Transport.Rd
 5ef7eeb790f9d686812407395a93dda9 *man/dispersion.Rd
-6384c0365cb1b46ea7230227dd9c8392 *man/electors.Rd
+19d83ee765d701464cbee7e9fc265abd *man/electors.Rd
 63859c49419644624290e47c206e1087 *man/getSummary-mclogit.Rd
-70fd8ab2f1d80443c6d03f65cec7bba2 *man/mblogit.Rd
-7371d905f3a1ee580c5165e026204362 *man/mclogit.Rd
+b3a28727c05a8cd22815f76f48dc2c98 *man/mblogit.Rd
+2f52366c045da9bf2712ca7a65c6d48d *man/mclogit.Rd
 e5473615d04e5ffc8f86cf2e5763b487 *man/mclogit.fit.Rd
 d5ede3a30233cb659a85d69bced029da *man/mclogit_control.Rd
-826973b468cfad32976855f741609ffb *man/simulate.Rd
+d77d5ecc5744e8a3f939e856f16724ae *man/predict.Rd
+1337900365bda2918af7aaa5cef8a210 *man/simulate.Rd

NAMESPACE

@@ -31,6 +31,7 @@ S3method(anova,mclogit)
 S3method(update,mclogit)
 
 S3method(print,mmclogit)
+S3method(vcov,mmclogit)
 S3method(summary,mmclogit)
 S3method(print,summary.mmclogit)
 
@@ -60,3 +61,6 @@ S3method(simulate,mblogit)
 
 S3method(simulate,mmclogit)
 S3method(simulate,mmblogit)
+
+S3method(predict,mmblogit)
+S3method(predict,mmclogit)

R/blockMatrices.R

@@ -84,9 +84,9 @@ bMatTrns <- function(x){
     m <- nrow(x)
     n <- ncol(x)
     res <- blockMatrix(nrow=n,ncol=m)
-    for(i in 1:m)
-        for(j in 1:n){
-            res[[i,j]] <- t(res[[i,j]])
+    for(i in 1:n)
+        for(j in 1:m){
+            res[[i,j]] <- t(x[[j,i]])
         }
     res
 }

R/mblogit.R

@@ -77,8 +77,8 @@
 #'    \url{https://doi.org/10.1080/01621459.1993.10594284}
 #'
 #' 
-#' @aliases predict.mblogit print.mblogit summary.mblogit print.summary.mblogit fitted.mblogit
-#' predict.mblogit, weights.mblogit
+#' @aliases print.mblogit summary.mblogit print.summary.mblogit fitted.mblogit
+#' weights.mblogit
 #' print.mmblogit summary.mmblogit print.summary.mmblogit
 mblogit <- function(formula,
                     data=parent.frame(),
@@ -174,6 +174,13 @@ mblogit <- function(formula,
             Y[cbind(i,j)] <- prior.weights
             w <- rowSums(Y)
             Y <- Y/w
+            if(any(w==0)){
+                Y[w==0,] <- 0
+                N <- sum(weights[w>0])
+                warning(sprintf("ignoring %d observerations with counts that sum to zero",
+                                sum(w==0)),
+                        call. = FALSE, immediate. = TRUE)
+            }
             Y <- as.vector(t(Y))
             weights <- rep(w,each=m)
             D <- diag(m)[,-1, drop=FALSE]
@@ -281,7 +288,8 @@ mblogit <- function(formula,
     }
     fit <- c(fit,list(call = call, formula = formula,
                       terms = mt,
-                      random = NULL,
+                      random = random,
+                      groups = groups,
                       data = data,
                       contrasts = contrasts,
                       xlevels = xlevels,
@@ -294,8 +302,9 @@ mblogit <- function(formula,
                       response.type=response.type,
                       from.table=from.table))
 
-    if(length(random))
+    if(length(random)){
         class(fit) <- c("mmblogit","mblogit","mmclogit","mclogit","lm")
+    }
     else
         class(fit) <- c("mblogit","mclogit","lm")
     fit
@@ -702,3 +711,159 @@ sample_factor <- function(probs, nsim =1, seed = NULL, ...){
     attr(yy,"seed") <- RNGstate
     return(yy)
 }
+
+lenuniq <- function(x) length(unique(x))
+
+predict.mmblogit <- function(object, newdata=NULL,type=c("link","response"),se.fit=FALSE,
+                             conditional=TRUE, ...){
+
+    type <- match.arg(type)
+    rhs <- object$formula[-2]
+    random <- object$random  
+    if(missing(newdata)){
+        mf <- object$model
+        na.act <- object$na.action
+    }
+    else{
+        vars <- unique(c(all.vars(rhs),all.vars(object$call$random),all.vars(object$call$weights)))
+        fo <- paste("~",paste(vars,collapse=" + "))
+        fo <- as.formula(fo,env=parent.frame())
+        mf <- model.frame(fo,data=newdata,na.action=na.exclude)
+        na.act <- attr(mf,"na.action")
+    }
+    X <- model.matrix(rhs,mf,
+                      contrasts.arg=object$contrasts,
+                      xlev=object$xlevels
+                      )
+    D <- object$D
+    XD <- X%x%D
+    eta <- c(XD %*% coef(object))
+
+    if(object$method=="PQL" && conditional){
+
+        rf <- random$formula
+        rt <- terms(rf)
+        groups <- random$groups
+        all.groups <- object$groups
+
+        Z <- model.matrix(rt,mf,
+                      contrasts.arg=object$contrasts,
+                      xlev=object$xlevels
+                      )
+        ZD <- Z%x%D
+
+        colnames(ZD) <- paste0(rep(colnames(D),ncol(Z)),
+                               "~",
+                               rep(colnames(Z),each=ncol(D)))
+        colnames(ZD) <- gsub("(Intercept)","1",colnames(ZD),fixed=TRUE)
+
+        groups <- mf[groups]
+        groups <- lapply(groups,as.integer)
+        nlev <- length(groups)
+        if(nlev > 1){
+            for(i in 2:nlev){
+                mm <- attr(all.groups[[i]],"unique")
+                mmm <- cumprod(mm)
+                groups[[i]] <- mmm[i]*groups[[i-1]]+groups[[i]]
+            }
+        }
+
+        groups <- lapply(groups,rep,each=nrow(D))
+
+        ZD <- Map(mkZ2,
+                  all.groups=all.groups,
+                  groups=groups,
+                  rX=list(ZD))
+        ZD <- blockMatrix(ZD)
+
+        random.effects <- object$random.effects
+        for(k in 1:nlev)
+            eta <- eta +  as.vector(ZD[[k]]%*%random.effects[[k]])
+    }
+
+    rspmat <- function(x){
+        y <- t(matrix(x,nrow=nrow(D)))
+        colnames(y) <- rownames(D)
+        y
+    }
+    eta <- rspmat(eta)
+
+    nvar <- ncol(X)
+    nobs <- nrow(X)
+
+    if(se.fit || type=="response"){
+        exp.eta <- exp(eta)
+        sum.exp.eta <- rowSums(exp.eta)
+        p <- exp.eta/sum.exp.eta
+    }
+    if(se.fit){
+        ncat <- ncol(p)
+        W <- Matrix(0,nrow=nobs*ncat,ncol=nobs)
+        i <- seq.int(ncat*nobs)
+        j <- rep(1:nobs,each=ncat)
+        pv <- as.vector(t(p))
+        W[cbind(i,j)] <- pv
+        W <- Diagonal(x=pv)-tcrossprod(W)
+        WX <- W%*%XD
+        if(object$method=="PQL"){
+            WZ <- bMatProd(W,ZD)
+            H <- object$info.fixed.random
+            K <- solve(H)
+        }
+    }
+
+    if(type=="response") {
+        if(se.fit){
+            if(object$method=="PQL" && conditional){
+                WXZ <- structure(cbind(blockMatrix(WX),WZ),class="blockMatrix")
+                var.p <- bMatProd(WXZ,K)
+                var.p <- Map(`*`,WXZ,var.p)
+                var.p <- lapply(var.p,rowSums)
+                var.p <- Reduce(`+`,var.p)
+            }
+            else {
+                vcov.coef <- vcov(object)
+                var.p <- rowSums(WX*(WX%*%vcov.coef))
+            }
+            se.p <- sqrt(var.p)
+            se.p <- rspmat(se.p)
+            if(is.null(na.act))
+                list(fit=p,se.fit=se.p) 
+            else
+                list(fit=napredict(na.act,p),
+                     se.fit=napredict(na.act,se.p))
+        }
+        else{
+            if(is.null(na.act)) p
+            else napredict(na.act,p)
+        }
+    }
+    else {
+        eta <- eta[,-1,drop=FALSE]
+        if(se.fit){
+            if(object$method=="PQL" && conditional){
+                XZ <- structure(cbind(blockMatrix(XD),ZD),class="blockMatrix")
+                var.eta <- bMatProd(XZ,K)
+                var.eta <- Map(`*`,XZ,var.eta)
+                var.eta <- lapply(var.eta,rowSums)
+                var.eta <- Reduce(`+`,var.eta)
+            }
+            else {
+                vcov.coef <- vcov(object)
+                var.eta <- rowSums(XD*(XD%*%vcov.coef))
+            }
+            se.eta <- sqrt(var.eta)
+            se.eta <- rspmat(se.eta)
+            se.eta <- se.eta[,-1,drop=FALSE]
+            if(is.null(na.act))
+                list(fit=eta,se.fit=se.eta) 
+            else
+                list(fit=napredict(na.act,eta),
+                     se.fit=napredict(na.act,se.eta))
+        }
+        else {
+            if(is.null(na.act)) eta
+            else napredict(na.act,eta)
+        }
+    }
+}

R/mclogit-dispersion.R

@@ -24,7 +24,7 @@ mclogit.dispersion <- function(y,w,s,pi,coef,method){
 
 update_mclogit_dispersion <- function(object,dispersion){
 
-    if(!isFALSE(dispersion)){
+    if(!missing(dispersion)){
         if(is.numeric(dispersion))
             phi <- dispersion
         else {

R/mclogit-fit.R

@@ -114,7 +114,6 @@ mclogit.fit <- function(
         phi = phi,
         linear.predictors = eta,
         working.residuals = (y-pi)/pi,
-        working.weights = w,
         response.residuals = y-pi,
         df.residual = resid.df,
         model.df = model.df,