[Depmix-commits] r225 - trunk/R

[noreply at r-forge.r-project.org]

Author: maarten
Date: 2008-08-26 23:39:56 +0200 (Tue, 26 Aug 2008)
New Revision: 225

Modified:
   trunk/R/responseGLM.R
   trunk/R/responseGLMMULTINOM.R
Log:
- fixed problem with n>1 in multinomial fitting and simulation
- fixed related (potential) bug in binomial

Modified: trunk/R/responseGLM.R
===================================================================
--- trunk/R/responseGLM.R	2008-08-26 13:16:23 UTC (rev 224)
+++ trunk/R/responseGLM.R	2008-08-26 21:39:56 UTC (rev 225)
@@ -35,33 +35,26 @@
 		if(family$family=="binomial") {
 			# FIX ME
 			y <- model.response(mf)
-			switch(is(y)[1],
-				factor = {
-					y <- as.matrix(as.numeric(as.numeric(y)==1))
-					#n <- matrix(1,nrow=nrow(y))
-				},
-				matrix = {
-					if(ncol(y) == 2) {
-						#n <- as.matrix(rowSums(y))
-						#y <- as.matrix(y[,1])
-					} else {
-						stop("model response not valid for binomial model")
-					}
-				},
-				numeric = {
-					if(sum(y %in% c(0,1)) != length(y)) stop("model response not valid for binomial model")
-					#n <- matrix(1,nrow=length(y))
+			if(NCOL(y) == 1) {
+        if(is.factor(y)) y <- as.matrix(as.numeric(as.numeric(y)==1)) else {
+          if(!is.numeric(y)) stop("model response not valid for binomial model")
+          if(sum(y %in% c(0,1)) != length(y)) stop("model response not valid for binomial model")
 					y <- as.matrix(y)
-				},
-				stop("model response not valid for binomial model")
-				# assume 1 success, rest not
-				#y <- as.numeric(as.numeric(y)==1)
-			)
+		   }
+			} else {
+			 if(ncol(y) != 2) {
+			   stop("model response not valid for binomial model")
+			 }
+			}
 		}
 		if(family$family=="multinomial") {
 			y <- model.response(mf)
-			# FIX THIS: comment out when it produces errors when the response is a matrix, see multinom help page
- 			if(is.factor(y)) y <- model.matrix(~y-1) else if(is.numeric(y)) y <- model.matrix(~factor(y)-1)
+			if(NCOL(y) == 1) {
+        if(is.factor(y)) y <- model.matrix(~y-1) else {
+          if(!is.numeric(y)) stop("model response not valid for binomial model")
+					y <- model.matrix(~factor(y)-1)
+		   }
+			}
 			parameters$coefficients <- matrix(0,ncol=ncol(y),nrow=ncol(x))
 			if(is.null(fixed)) {
 				fixed <- parameters$coefficients
Modified: trunk/R/responseGLMMULTINOM.R
===================================================================
--- trunk/R/responseGLMMULTINOM.R	2008-08-26 13:16:23 UTC (rev 224)
+++ trunk/R/responseGLMMULTINOM.R	2008-08-26 21:39:56 UTC (rev 225)
@@ -61,12 +61,14 @@
 		if(missing(times)) {
 			# draw all times in one go
 			pr <- predict(object)
+			n <- rowSums(object at y)
 		} else {
 			pr <- predict(object)[times,]
+			n <- rowSums(object at y)[times]
 			if(length(times)==1) pr <- matrix(pr,ncol=length(pr))
 		}
 		nt <- nrow(pr)
-		sims <- array(apply(pr,1,rmultinom,n=nsim,size=1),dim=c(ncol(pr),nsim,nt))
+		sims <- array(apply(pr,1,rmultinom,n=nsim,size=n),dim=c(ncol(pr),nsim,nt))
 		sims <- matrix(aperm(sims,c(3,2,1)),nrow=nsim*nt,ncol=ncol(pr))
 		#response <- t(apply(sims,c(2,3), function(x) which(x==1)))
 		return(sims)