[Depmix-commits] r317 - trunk/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 19 15:30:43 CET 2010
Author: ingmarvisser
Date: 2010-01-19 15:30:42 +0100 (Tue, 19 Jan 2010)
New Revision: 317
Modified:
trunk/R/multinomial.R
trunk/R/responseGLM.R
trunk/R/responseGLMMULTINOM.R
trunk/R/responseNORM.R
trunk/R/transInit.R
Log:
Added multinomial response function with identity link (no covariates allowed).
Modified: trunk/R/multinomial.R
===================================================================
--- trunk/R/multinomial.R 2010-01-19 12:37:01 UTC (rev 316)
+++ trunk/R/multinomial.R 2010-01-19 14:30:42 UTC (rev 317)
@@ -20,11 +20,11 @@
if (linktemp %in% okLinks) {
if(linktemp == "mlogit") stats <- mlogit() else stats <- make.link(linktemp)
} else {
- if (is.character(link)) {
+ if(is.character(link)) {
stats <- make.link(link)
linktemp <- link
} else {
- if (inherits(link, "link-glm")) {
+ if(inherits(link, "link-glm")) {
stats <- link
if (!is.null(stats$name))
linktemp <- stats$name
@@ -42,13 +42,13 @@
validmu <- function(mu) {
all(mu > 0) && all(mu < 1)
}
- dev.resids <- function(y,mu,wt) {
-
+
+ dev.resids <- function(y,mu,wt) {
}
initialize <- expression()
structure(list(family = "multinomial", link = linktemp, linkfun = stats$linkfun,
linkinv = stats$linkinv, variance = variance, dev.resids = dev.resids,
mu.eta = stats$mu.eta, initialize = initialize, validmu = validmu, valideta = stats$valideta, base=base),
- class = "family")
+ class = "family")
}
Modified: trunk/R/responseGLM.R
===================================================================
--- trunk/R/responseGLM.R 2010-01-19 12:37:01 UTC (rev 316)
+++ trunk/R/responseGLM.R 2010-01-19 14:30:42 UTC (rev 317)
@@ -50,15 +50,25 @@
if(family$family=="multinomial") {
y <- model.response(mf)
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")
+ if(is.factor(y)) {
+ y <- model.matrix(~y-1)
+ } else {
+ if(!is.numeric(y)) stop("model response not valid for multinomial model")
y <- model.matrix(~factor(y)-1)
}
}
- parameters$coefficients <- matrix(0,ncol=ncol(y),nrow=ncol(x))
- if(is.null(fixed)) {
+ if(family$link=="mlogit") {
+ parameters$coefficients <- matrix(0,ncol=ncol(y),nrow=ncol(x))
+ if(is.null(fixed)) {
+ fixed <- parameters$coefficients
+ fixed[,family$base] <- 1
+ fixed <- c(as.logical(t(fixed)))
+ }
+ }
+ 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),nrow=ncol(x))
fixed <- parameters$coefficients
- fixed[,family$base] <- 1
fixed <- c(as.logical(t(fixed)))
}
}
@@ -76,16 +86,16 @@
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)])
- }
+ 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,
@@ -106,7 +116,7 @@
print(object at formula)
cat("Coefficients: \n")
print(object at parameters$coefficients)
- if(object at family$family=="multinomial") {
+ if(object at family$family=="multinomial"&object at family$link!="identity") {
# also print probabilities at covariate values of zero
cat("Probalities at zero values of the covariates.\n")
if(!(is.null(dim(object at parameters$coefficients)))) {
Modified: trunk/R/responseGLMMULTINOM.R
===================================================================
--- trunk/R/responseGLMMULTINOM.R 2010-01-19 12:37:01 UTC (rev 316)
+++ trunk/R/responseGLMMULTINOM.R 2010-01-19 14:30:42 UTC (rev 317)
@@ -4,34 +4,42 @@
setMethod("fit","MULTINOMresponse",
function(object,w) {
if(missing(w)) w <- NULL
- pars <- object at parameters
- base <- object at family$base # delete me
- y <- object at y
- x <- object at x
- #if(is.null(w)) w <- rep(1,nrow(y))
- # mask is an nx*ny matrix (x are inputs, y are output levels)
- mask <- matrix(1,nrow=nrow(pars$coefficients),ncol=ncol(pars$coefficients))
- mask[,base] <- 0 # fix base category coefficients to 0
- mask <- rbind(0,mask) # fix automatic "bias" nodes to 0
- Wts <- mask
- Wts[-1,] <- pars$coefficients # set starting weights
- if(!is.null(w)) {
- if(NCOL(y) < 3) {
- fit <- nnet.default(x,y,weights=w,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+ if(object at family$link=="mlogit") {
+ pars <- object at parameters
+ base <- object at family$base # delete me
+ y <- object at y
+ x <- object at x
+ #if(is.null(w)) w <- rep(1,nrow(y))
+ # mask is an nx*ny matrix (x are inputs, y are output levels)
+ mask <- matrix(1,nrow=nrow(pars$coefficients),ncol=ncol(pars$coefficients))
+ mask[,base] <- 0 # fix base category coefficients to 0
+ mask <- rbind(0,mask) # fix automatic "bias" nodes to 0
+ Wts <- mask
+ Wts[-1,] <- pars$coefficients # set starting weights
+ if(!is.null(w)) {
+ if(NCOL(y) < 3) {
+ fit <- nnet.default(x,y,weights=w,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+ } else {
+ fit <- nnet.default(x,y,weights=w,size=0,softmax=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+ }
} else {
- fit <- nnet.default(x,y,weights=w,size=0,softmax=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+ if(NCOL(y) < 3) {
+ fit <- nnet.default(x,y,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+ } else {
+ fit <- nnet.default(x,y,size=0,softmax=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+ }
}
- } else {
- if(NCOL(y) < 3) {
- fit <- nnet.default(x,y,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
- } else {
- fit <- nnet.default(x,y,size=0,softmax=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
- }
+ # this is necessary because setpars wants coefficients in column major order
+ pars$coefficients <- t(matrix(fit$wts,ncol=ncol(pars$coefficients),nrow=nrow(pars$coefficients)+1)[-1,])
+ # parameters are set correctly now
+ object <- setpars(object,unlist(pars))
}
- # this is necessary because setpars wants coefficients in column major order
- pars$coefficients <- t(matrix(fit$wts,ncol=ncol(pars$coefficients),nrow=nrow(pars$coefficients)+1)[-1,])
- # parameters are set correctly now
- object <- setpars(object,unlist(pars))
+ if(object at family$link=="identity") {
+ if(is.null(w)) w <- rep(1,nrow(object at y))
+ sw <- sum(w)
+ pars <- c(apply(object at y,2,function(x){sum(x*w)/sw}))
+ object <- setpars(object,pars)
+ }
object
}
)
Modified: trunk/R/responseNORM.R
===================================================================
--- trunk/R/responseNORM.R 2010-01-19 12:37:01 UTC (rev 316)
+++ trunk/R/responseNORM.R 2010-01-19 14:30:42 UTC (rev 317)
@@ -7,19 +7,19 @@
setMethod("fit","NORMresponse",
function(object,w) {
- if(missing(w)) w <- NULL
+ if(missing(w)) w <- NULL
pars <- object at parameters
if(!is.null(w)) {
- fit <- lm.wfit(x=object at x,y=object at y,w=w)
- } else {
- fit <- lm.fit(x=object at x,y=object at y)
- }
+ fit <- lm.wfit(x=object at x,y=object at y,w=w)
+ } else {
+ fit <- lm.fit(x=object at x,y=object at y)
+ }
pars$coefficients <- fit$coefficients
if(!is.null(w)) {
- pars$sd <- sqrt(sum(w*fit$residuals^2/sum(w)))
- } else {
- pars$sd <- sd(fit$residuals)
- }
+ pars$sd <- sqrt(sum(w*fit$residuals^2/sum(w)))
+ } else {
+ pars$sd <- sd(fit$residuals)
+ }
object <- setpars(object,unlist(pars))
object
}
Modified: trunk/R/transInit.R
===================================================================
--- trunk/R/transInit.R 2010-01-19 12:37:01 UTC (rev 316)
+++ trunk/R/transInit.R 2010-01-19 14:30:42 UTC (rev 317)
@@ -78,8 +78,8 @@
if(object at family$link=="identity") object at family$linkinv(object at x%*%object at parameters$coefficients)
else object at family$linkinv(object at x%*%object at parameters$coefficients,base=object at family$base)
} else {
- if(!(is.matrix(newx))) stop("'newx' must be matrix in predict (GLMresponse)")
- if(!(ncol(newx)==nrow(object at parameters$coefficients))) stop("Incorrect dimension of 'newx' in predict (GLMresponse)")
+ if(!(is.matrix(newx))) stop("'newx' must be matrix in predict(transInit)")
+ if(!(ncol(newx)==nrow(object at parameters$coefficients))) stop("Incorrect dimension of 'newx' in predict(transInit)")
if(object at family$link=="identity") object at family$linkinv(newx%*%object at parameters$coefficients)
else object at family$linkinv(newx%*%object at parameters$coefficients,base=object at family$base)
}
More information about the depmix-commits
mailing list