[Depmix-commits] r552 - pkg/depmixS4/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jul 30 14:33:44 CEST 2012
Author: ingmarvisser
Date: 2012-07-30 14:33:44 +0200 (Mon, 30 Jul 2012)
New Revision: 552
Modified:
pkg/depmixS4/R/responseGLM.R
Log:
- changed parameters$coefficients from matrix to vector in multinomial identity models
- added names to those
- changed getpars such that it preserves names of parameters
- changed setpars such that it preserves names of parameters
Modified: pkg/depmixS4/R/responseGLM.R
===================================================================
--- pkg/depmixS4/R/responseGLM.R 2012-07-26 13:46:55 UTC (rev 551)
+++ pkg/depmixS4/R/responseGLM.R 2012-07-30 12:33:44 UTC (rev 552)
@@ -14,6 +14,7 @@
contains="response"
)
+
setMethod("GLMresponse",
signature(formula="formula"),
function(formula,data=NULL,family=gaussian(),pstart=NULL,fixed=NULL,prob=TRUE,na.action="na.pass", ...) {
@@ -32,6 +33,7 @@
parameters <- list()
constr <- NULL
parameters$coefficients <- vector("numeric",length=ncol(x))
+ names(parameters$coefficients) <- attr(x,"dimnames")[[2]]
if(family$family=="gaussian") {
parameters$sd <- 1
constr <- list(
@@ -58,11 +60,13 @@
y <- model.response(mf)
if(NCOL(y) == 1) {
if(is.factor(y)) {
+ namesy <- levels(y)
mf <- model.frame(~y-1,na.action=na.action)
y <- model.matrix(attr(mf, "terms"),mf)
#y <- model.matrix(~y-1,na.action=na.action)
} else {
if(!is.numeric(y)) stop("model response not valid for multinomial model")
+ namesy <- levels(factor(y))
mf <- model.frame(~factor(y)-1,na.action=na.action)
y <- model.matrix(attr(mf, "terms"),mf)
#y <- model.matrix(~factor(y)-1,na.action=na.action)
@@ -75,11 +79,16 @@
fixed <- parameters$coefficients
fixed[,family$base] <- 1
fixed <- c(as.logical(t(fixed)))
- }
+ }
+ if(is.null(namesy)) namesy <- 1:ncol(y)
+ colnames(parameters$coefficients) <- namesy
+ rownames(parameters$coefficients) <- attr(x,"dimnames")[[2]]
}
if(family$link=="identity") {
if(ncol(x)>1) stop("covariates not allowed in multinomial model with identity link")
- parameters$coefficients <- matrix(1/ncol(y),ncol=ncol(y),1)
+ ncy <- ncol(y)
+ parameters$coefficients <- rep(1/ncy,ncy)
+ names(parameters$coefficients) <- paste("pr",1:ncol(y),sep="")
fixed <- rep(0,ncol(y))
fixed <- c(as.logical(t(fixed)))
constr <- list(
@@ -88,49 +97,50 @@
linlow = 1,
parup = rep(1,ncol(y)),
parlow = rep(0,ncol(y))
- )
- }
+ )
+ }
}
npar <- length(unlist(parameters))
if(is.null(fixed)) fixed <- as.logical(rep(0,npar))
if(!is.null(pstart)) {
- if(length(pstart)!=npar) stop("length of 'pstart' must be",npar)
- if(family$family=="multinomial") {
- if(family$link=="identity") {
- parameters$coefficients[1,] <- pstart[1:ncol(parameters$coefficients)]/sum(pstart[1:ncol(parameters$coefficients)])
-
+ if(length(pstart)!=npar) stop("length of 'pstart' must be",npar)
+ if(family$family=="multinomial") {
+ if(family$link=="identity") {
+ parameters$coefficients <- pstart[1:length(parameters$coefficients)]/sum(pstart[1:length(parameters$coefficients)])
+ } else {
+ if(prob) parameters$coefficients[1,] <- family$linkfun(pstart[1:ncol(parameters$coefficients)],base=family$base)
+ 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),]
} else {
- if(prob) parameters$coefficients[1,] <- family$linkfun(pstart[1:ncol(parameters$coefficients)],base=family$base)
- else parameters$coefficients[1,] <- pstart[1:ncol(parameters$coefficients)]
+ # if(prob) parameters$coefficients <- family$linkfun(as.numeric(pstart[1:length(parameters$coefficients)]))
+ if(family$family=="binomial") {
+ if(prob) parameters$coefficients[1] <- family$linkfun(pstart[1])
+ else parameters$coefficients[1] <- pstart[1]
+ if(ncol(x)>1) parameters$coefficients[2:ncol(x)] <- pstart[2:ncol(x)]
+ } else {
+ parameters$coefficients[1:ncol(x)] <- pstart[1:ncol(x)]
+ }
+ if(length(unlist(parameters))>length(parameters$coefficients)) {
+ if(family$family=="gaussian") parameters$sd <- as.numeric(pstart[(length(parameters$coefficients)+1)])
+ }
}
- 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 <- family$linkfun(as.numeric(pstart[1:length(parameters$coefficients)]))
- if(family$family=="binomial") {
- if(prob) parameters$coefficients[1] <- family$linkfun(pstart[1])
- else parameters$coefficients[1] <- pstart[1]
- if(ncol(x)>1) parameters$coefficients[2:ncol(x)] <- pstart[2:ncol(x)]
- } else {
- parameters$coefficients[1:ncol(x)] <- pstart[1:ncol(x)]
- }
- if(length(unlist(parameters))>length(parameters$coefficients)) {
- if(family$family=="gaussian") parameters$sd <- as.numeric(pstart[(length(parameters$coefficients)+1)])
- }
- }
}
mod <- switch(family$family,
- gaussian = new("NORMresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
- binomial = new("BINOMresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
- multinomial = new("MULTINOMresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
- poisson = new("POISSONresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
- Gamma = new("GAMMAresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
- new("GLMresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr)
+ gaussian = new("NORMresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
+ binomial = new("BINOMresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
+ multinomial = new("MULTINOMresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
+ poisson = new("POISSONresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
+ Gamma = new("GAMMAresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr),
+ new("GLMresponse",formula=formula,family=family,parameters=parameters,fixed=fixed,x=x,y=y,npar=npar,constr=constr)
)
mod
- }
+ }
)
+
+
setMethod("show","GLMresponse",
function(object) {
cat("Model of type ", object at family$family," (", object at family$link, "), formula: ", sep="")
@@ -169,7 +179,7 @@
npar <- npar(object)
if(length(values)!=npar) stop("length of 'values' must be",npar)
# determine whether parameters or fixed constraints are being set
- nms <- names(object at parameters$coefficients)
+ nms <- attributes(object at parameters$coefficients)
if(length(values) == 0) return(object) # nothing to set;
switch(which,
"pars"= {
@@ -188,31 +198,33 @@
"fixed" = {
object at fixed <- as.logical(values)
}
- )
- names(object at parameters$coefficients) <- nms
+ )
+ attributes(object at parameters$coefficients) <- nms
return(object)
}
)
+
setMethod("getpars","GLMresponse",
- function(object,which="pars",...) {
- switch(which,
- "pars" = {
- parameters <- numeric()
- if(object at family$family=="multinomial") {
- # coefficient is usually a matrix here
- parameters <- c(t(object at parameters$coefficients)) # Why transpose?
- } else {
- parameters <- unlist(object at parameters)
- }
- pars <- parameters
- },
- "fixed" = {
- pars <- object at fixed
- }
- )
- return(pars)
- }
+ function(object,which="pars",...) {
+ switch(which,
+ "pars" = {
+ parameters <- numeric()
+ if(object at family$family=="multinomial"&object at family$link=="mlogit") {
+ # coefficient is usually a matrix here
+ parameters <- c(t(object at parameters$coefficients)) # Why transpose?
+ } else {
+ parameters <- object at parameters$coefficients
+ if(object at family$family=="gaussian") parameters <- c(parameters,object at parameters$sd)
+ }
+ pars <- parameters
+ },
+ "fixed" = {
+ pars <- object at fixed
+ }
+ )
+ return(pars)
+ }
)
# methods: fit, logDens, predict
More information about the depmix-commits
mailing list