[Depmix-commits] r381 - trunk/R

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

Author: maarten
Date: 2010-03-04 18:10:32 +0100 (Thu, 04 Mar 2010)
New Revision: 381

Modified:
   trunk/R/transInit.R
Log:
- corrected starting values for transInit models with identity links
- transInit method no longer needs a data.frame
- TODO: AIC and BIC are not computed correctly for transInit models with identity link!

Modified: trunk/R/transInit.R
===================================================================
--- trunk/R/transInit.R	2010-03-04 16:05:35 UTC (rev 380)
+++ trunk/R/transInit.R	2010-03-04 17:10:32 UTC (rev 381)
@@ -11,41 +11,55 @@
 setMethod("transInit",
 	signature(formula="formula"),
 	function(formula,nstates,data=NULL,family=multinomial(),pstart=NULL,fixed=NULL,prob=TRUE, ...) {
-		call <- match.call()
-		mf <- match.call(expand.dots = FALSE)
-		m <- match(c("formula", "data"), names(mf), 0)
-		mf <- mf[c(1, m)]
-		mf$drop.unused.levels <- TRUE
-		mf[[1]] <- as.name("model.frame")
-		mf <- eval(mf, parent.frame())
-		x <- model.matrix(attr(mf, "terms"),mf)
+		call <- match.call()
+		if(formula==formula(~1) & is.null(data)) {
+			x <- matrix(1,ncol=1)
+		} else {
+			mf <- match.call(expand.dots = FALSE)
+			m <- match(c("formula", "data"), names(mf), 0)
+			mf <- mf[c(1, m)]
+			mf$drop.unused.levels <- TRUE
+			mf[[1]] <- as.name("model.frame")
+			mf <- eval(mf, parent.frame())
+			x <- model.matrix(attr(mf, "terms"),mf)
+		}
 		y <- matrix(1,ncol=1) # y is not needed in the transition and init models
 		parameters <- list()
 		if(is.null(nstates)) stop("'nstates' must be provided in call to trinModel")
-		if(family$family=="multinomial") {
-			parameters$coefficients <- matrix(0,ncol=nstates,nrow=ncol(x))
-			if(is.null(fixed)) {
-				fixed <- parameters$coefficients
-				fixed[,family$base] <- 1 
-				fixed <- c(as.logical(t(fixed)))
+		if(family$family=="multinomial") {
+			if(family$link=="identity") {
+				parameters$coefficients <- t(apply(matrix(1,ncol=nstates,nrow=ncol(x)),1,function(x) x/sum(x)))
+				if(is.null(fixed)) {
+					fixed <- matrix(0,nrow=nrow(parameters$coefficients),ncol=ncol(parameters$coefficients))
+					fixed <- c(as.logical(t(fixed)))
+				}
+			} else {
+				parameters$coefficients <- matrix(0,ncol=nstates,nrow=ncol(x))
+				if(is.null(fixed)) {
+					fixed <- parameters$coefficients
+					fixed[,family$base] <- 1 
+					fixed <- c(as.logical(t(fixed)))
+				}
 			}
 		}
 		npar <- length(unlist(parameters))
 		if(is.null(fixed)) fixed <- rep(0,npar)
 		if(!is.null(pstart)) {
 			if(length(pstart)!=npar) stop("length of 'pstart' must be ",npar)
-			if(family$family=="multinomial") {
-				if(prob) {
-					if(family$link=="identity") {
-						parameters$coefficients[1,] <- family$linkfun(pstart[1:ncol(parameters$coefficients)])
+			if(family$family=="multinomial") {
+				if(family$link=="identity") {
+					parameters$coefficients[1,] <- pstart[1:ncol(parameters$coefficients)]
+					pstart <- matrix(pstart,ncol(x),byrow=TRUE)
+					if(ncol(x)>1) parameters$coefficients[2:ncol(x),] <- pstart[2:ncol(x),]
+				} else {
+					if(prob) {
+						parameters$coefficients[1,] <- family$linkfun(pstart[1:ncol(parameters$coefficients)],base=family$base)
 					} else {
-						parameters$coefficients[1,] <- family$linkfun(pstart[1:ncol(parameters$coefficients)],base=family$base)
+						parameters$coefficients[1,] <- pstart[1:ncol(parameters$coefficients)]
 					}
-				} else {
-					parameters$coefficients[1,] <- pstart[1:ncol(parameters$coefficients)]
+					pstart <- matrix(pstart,ncol(x),byrow=TRUE)
+					if(ncol(x)>1) parameters$coefficients[2:ncol(x),] <- pstart[2:ncol(x),]
 				}
-				pstart <- matrix(pstart,ncol(x),byrow=TRUE)
-				if(ncol(x)>1) parameters$coefficients[2:ncol(x),] <- pstart[2:ncol(x),]
 			} else {
 				if(family$link=="identity") parameters$coefficients <- family$linkfun(pstart[1:length(parameters$coefficients)])
 				else parameters$coefficients <- family$linkfun(pstart[1:length(parameters$coefficients)],base=family$base)
@@ -176,3 +190,13 @@
 		}
 	}
 )
+
+setMethod("getdf","transInit",
+	function(object) {
+		npar <- sum(!object at fixed)
+		if(object at family$linkfun == "identity") {
+			fx <- t(matrix(object at fixed,ncol=nrow(object at parameters$coefficients),nrow=ncol(object at parameters$coefficients)))
+			npar <- npar - rowSums(fx)
+		}
+	}
+)