[Depmix-commits] r171 - trunk/R

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

Author: ingmarvisser
Date: 2008-06-17 16:11:27 +0200 (Tue, 17 Jun 2008)
New Revision: 171

Modified:
   trunk/R/responseGLM.R
   trunk/R/responseGLMBINOM.R
Log:
Fixed a bug in dens function for binomial models

Modified: trunk/R/responseGLM.R
===================================================================
--- trunk/R/responseGLM.R	2008-06-17 12:52:04 UTC (rev 170)
+++ trunk/R/responseGLM.R	2008-06-17 14:11:27 UTC (rev 171)
@@ -34,32 +34,32 @@
 		}
 		if(family$family=="binomial") {
 			# FIX ME
-		  y <- model.response(mf)
+			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))
-    		  y <- as.matrix(y)
-    		},
-    		stop("model response not valid for binomial model")
-    			 # assume 1 success, rest not
-    			 #y <- as.numeric(as.numeric(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))
+					y <- as.matrix(y)
+				},
+				stop("model response not valid for binomial model")
+				# assume 1 success, rest not
+				#y <- as.numeric(as.numeric(y)==1)
+			)
+		}
 		if(family$family=="multinomial") {
-      y <- model.response(mf)
+			y <- model.response(mf)
 			if(is.factor(y)) y <- model.matrix(~y-1) else if(is.numeric(y)) y <- model.matrix(~factor(y)-1)
 			parameters$coefficients <- matrix(0,ncol=ncol(y),nrow=ncol(x))
 			if(is.null(fixed)) {
@@ -105,6 +105,16 @@
 		print(object at formula)
 		cat("Coefficients: \n")
 		print(object at parameters$coefficients)
+		if(object at family$family=="multinomial") {
+			# also print probabilities at covariate values of zero
+			cat("Probalities at zero values of the covariates.\n")
+			cat(object at family$linkinv(object at parameters$coefficients[1,],base=object at family$base),"\n")
+		}
+		if(object at family$family=="binomial") {
+			# also print probabilities at covariate values of zero
+			cat("Probality at zero values of the covariates.","\n")
+			cat(object at family$linkinv(object at parameters$coefficients[1]),"\n")
+		}
 		if(object at family$family=="gaussian") {
 			cat("sd ",object at parameters$sd,"\n")
 		}	
Modified: trunk/R/responseGLMBINOM.R
===================================================================
--- trunk/R/responseGLMBINOM.R	2008-06-17 12:52:04 UTC (rev 170)
+++ trunk/R/responseGLMBINOM.R	2008-06-17 14:11:27 UTC (rev 171)
@@ -10,42 +10,42 @@
 # returns: matrix with log(p(y|x,parameters))
 setMethod("logDens","BINOMresponse",
 	function(object) {
-    if(NCOL(object at y) == 2) {
-      dbinom(x=object at y[,1],size=object at y[,2],prob=predict(object),log=TRUE)
-    } else {
-      if(!NCOL(object at y==1)) stop("not a valid response matrix for BINOMresponse")
-		  dbinom(x=object at y,prob=predict(object),log=TRUE)
-    }
+		if(NCOL(object at y) == 2) {
+			dbinom(x=object at y[,1],size=rowSums(object at y),prob=predict(object),log=TRUE)
+		} else {
+			if(!NCOL(object at y==1)) stop("not a valid response matrix for BINOMresponse")
+			dbinom(x=object at y,prob=predict(object),log=TRUE)
+		}
 	}
 )
 
 setMethod("dens","BINOMresponse",
 	function(object,log=FALSE) {
-    if(NCOL(object at y) == 2) {
-      dbinom(x=object at y[,1],size=object at y[,2],prob=predict(object),log=log)
-    } else {
-      if(!NCOL(object at y==1)) stop("not a valid response matrix for BINOMresponse")
-		  dbinom(x=object at y,prob=predict(object),log=log)
-    }
+		if(NCOL(object at y) == 2) {
+			dbinom(x=object at y[,1],size=rowSums(object at y),prob=predict(object),log=log)
+		} else {
+			if(!NCOL(object at y==1)) stop("not a valid response matrix for BINOMresponse")
+			dbinom(x=object at y,prob=predict(object),log=log)
+		}
 	}
 )
 
 setMethod("simulate",signature(object="BINOMresponse"),
-  function(object,nsim=1,seed=NULL,time) {
-    if(missing(time)) {
-      # draw in one go
-      pr <- predict(object)
-    } else {
-      pr <- predict(object)[time,]
-    }
-    nt <- nrow(pr)
-    if(NCOL(object at y) == 2) {
-      response <- rbinom(nt*nsim,size=object at y[,2],prob=pr)
-    } else {
-      response <- rbinom(nt*nsim,size=1,prob=pr)
-    }
-    #if(nsim > 1) response <- matrix(response,ncol=nsim)
-    response <- as.matrix(response)
-    return(response)
-  }
+	function(object,nsim=1,seed=NULL,time) {
+		if(missing(time)) {
+			# draw in one go
+			pr <- predict(object)
+		} else {
+			pr <- predict(object)[time,]
+		}
+		nt <- nrow(pr)
+		if(NCOL(object at y) == 2) {
+			response <- rbinom(nt*nsim,size=object at y[,2],prob=pr)
+		} else {
+			response <- rbinom(nt*nsim,size=1,prob=pr)
+		}
+		#if(nsim > 1) response <- matrix(response,ncol=nsim)
+		response <- as.matrix(response)
+		return(response)
+	}
 )