R/GlmResponse.R

# GLM Response Model
### A response model which can be fitted by glm.fit
setClass("GlmResponse",
  contains="ResponseModel",
  representation(
    family="ANY",
    #formula="formula",
    #data = "list",
    z = "matrix",
    combine = "function"
    #LearnModPredName = "character"
    #sigma="matrix"
  )
)

setMethod("getReplication",signature=(object="GlmResponse"),
  function(object,case,ntimes=NULL,...) {
    if(length(object@z) > 0) {
    if(is.null(ntimes)) {
      lt <- object@nTimes@cases
      bt <- object@nTimes@bt
      et <- object@nTimes@et
    } else {
      lt <- length(ntimes)
      et <- cumsum(ntimes)
      bt <- c(1,et[-lt]+1)
    }
    z <- object@z[bt[case]:et[case],]
    } else {
      z <- object@z
    }
    out <- callNextMethod()
    out$z <- z
    out
  }
)

setMethod("fit",signature(object="GlmResponse"),
	function(object,...) {
    if(object@nTimes@cases > 1 && !object@parStruct@replicate) {
      if(length(object@parStruct@id)>0) stop("Constraints in GlmResponse models are currently not implemented")
      for(case in 1:object@nTimes@cases) {
        repl <- getReplication(object,case=case,...)
        x <- repl$x
        y <- repl$y
        #z <- repl$z
        pars <- repl$parameters
        #if(length(z) > 0) {
        #  x <- object@combine(x,z)
        #}
        runm <- glm.fit(x=x,y=y,family=object@family,start=pars$coefficients)
        pars$coefficients <- runm$coefficients
        if(object@family$family=="gaussian") pars$sd <- sqrt(sum((object@y - predict(object))^2)/(length(object@y)-1))
        object@parStruct@parameters[[case]] <- pars
      }
    } else {
      # y = response
      pars <- object@parStruct@parameters
      #if(length(objecct@z)>0) {
      #  x <- object@combine(object@x,object@z)
      #} else {
      #  x <- object@x
      #}
      runm <- glm.fit(x=object@x,y=object@y,family=object@family,start=pars$coefficients)
      pars$coefficients <- runm$coefficients
      if(object@family$family=="gaussian") pars$sd <- sqrt(sum((object@y - predict(object))^2)/(length(object@y)-1))
      object@parStruct@parameters <- setPars(object,unlist(pars),rval="parameters",...)
    }
    #object <- setpars(object,unlist(pars))
    object <- runm(object,...)
    return(object)
	}
)

setMethod("predict","GlmResponse",
	function(object,type="link") {
    if(object@nTimes@cases > 1 & !object@parStruct@replicate) {
      mu <- vector()
      for(case in 1:object@nTimes@cases) {
        if(NCOL(object@x) > 1) mu <- rbind(mu,object@x[object@nTimes@bt[case]:object@nTimes@et[case],]%*%as.matrix(object@parStruct@parameters[[case]]$coefficients)) else mu <- c(mu,object@parStruct@parameters[[case]]$coefficients*object@x[object@nTimes@bt[case]:object@nTimes@et[case],])
      }
    } else {
    # y = response
      if(NCOL(object@x) > 1) mu <- object@x%*%as.matrix(object@parStruct@parameters$coefficients) else mu <- object@parStruct@parameters$coefficients*object@x
    }
    if(type=="link") return(mu) else {
      if(type=="response") {
        return(object@family$linkinv(mu))
      }
    }
	}
)
setMethod("logLik","GlmResponse",
  function(object) {
    nt <- NROW(object@y)
    LL <- switch(object@family$family,
      gaussian = {
        mu <- predict(object,type="response")
        sum(dnorm(x=object@y,mean=mu,sd=object@parStruct@parameters$sd,log=TRUE))
      },
      binomial = {
        p <- predict(object,type="response")
        sum(dbinom(x=object@y,size=1,prob=p,log=TRUE))
      },
      poisson = {
        lambda <- predict(object,type="response")
        sum(dpois(x=object@y,lambda=lambda,log=TRUE))
      },
      Gamma = {
        shape <- predict(object,type="response")
        sum(dgamma(x=object@y,shape=shape,log=TRUE))
      },
      stop("family",object@family$family,"not implemented (yet)")
    )
    attr(LL,"nobs") <- nt
    attr(LL,"df") <- length(getPars(object,which="free"))
    class(LL) <- "logLik"
    LL
  }
)

setMethod("simulate",signature(object="GlmResponse"),
	function(object,nsim=1,seed=NULL,times) {
    if(!is.null(seed)) {
      old.seed <- .Random.seed 
      set.seed(seed)
    }
    if(missing(times)) {
      pr <- predict(object,type=response)
    } else {
      pr <- predict(object,type=response)[times,]
    }
    nt <- nrow(pr)
    response <- switch(object@family$family,
      gaussian = {
        rnorm(nt*nsim,mean=pr,sd=object@parStruct@parameters$sd)
      },
      binomial = {
        if(NCOL(object@y) == 2) {
    			rbinom(nt*nsim,size=object@y[,2],prob=pr)
    		} else {
    			rbinom(nt*nsim,size=1,prob=pr)
    		}
      },
      poisson = {
        rpois(nt*nsim,lambda=pr)
      },
      Gamma = {
        rgamma(nt*nsim,shape=pr)
      },
      stop("family",object@family$family,"not implemented (yet)")
    )
		#if(nsim > 1) response <- matrix(response,ncol=nsim)
		#response <- as.matrix(response)
		
		object@y <- as.matrix(response)
		if(!missing(times)) {
      object@x <- object@x[rep(times,nsim),]
      ntim <- rep(0,length=object@nTimes@cases)
  		for(i in 1:length(ntim)) {
  		  ntim[i] <- sum(seq(object@nTimes@bt[i],object@nTimes@et[i]) %in% times)
      }
      warning("simulation with a times argument may result in wrong parStruct argument; please check parameters.")
      object@parStruct <- rep(object@parStruct,times=nsim)
    } else {
      object@x <- object@x[rep(1:nrow(object@x),nsim),]
      ntim <- object@nTimes@n
      object@parStruct <- rep(object@parStruct,times=nsim)
    }
    ntim <- rep(ntim,nsim)
    object@nTimes <- nTimes(ntim)
    if(!is.null(seed)) {
      set.seed(old.seed)
    }
		return(object)
	}
)

setMethod("lFr",signature(x="LearningModel",y="GlmResponse"),
  function(x,y,...,type="link") {
    if(length(y@z) > 0) {
      xx <- predict(x,type=type,...)
      y@x <- y@combine(xx,y@z,...)
      #assign(as.character(y@LearnModPredName),predict(x,type="link",...))
      #form <- y@formula
      #form[[2]] <- NULL
      #y@x <- model.matrix(object=form,data=y@data,...)
    } else {
      y@x <- predict(x,type=type,...)
    }
    #y@x <- eval(y@x.call,envir=list(parent.frame(),y@data)) #predict(x,type="link",...)
    y
  }
)

GlmResponse <- function(formula,family=gaussian(),data,covariate,combine=function(x,z) { cbind(x,z) },parameters=list(),ntimes=NULL,replicate=TRUE,fixed,base=NULL,
                        parStruct,subset) {
  mf <- match.call(expand.dots = FALSE)
  m <- match(c("formula", "data", "subset","base"), names(mf), 0)
  mf <- mf[c(1, m)]
  mf$drop.unused.levels <- TRUE
  mf[[1]] <- as.name("model.frame")
  #mf[[2]] <- eval(substitute(mf$formula, list(as.name(predname) = as.name("learningModelPrediction"))))
  mf <- eval(mf, parent.frame())
  
  #if(!missing(subset)) dat <- mcpl.prepare(formula=formula,data=data,subset=subset,base=base) else
  #  dat <- mcpl.prepare(formula=formula,data=data,base=base)
  y <- dat$y
  x <- dat$x
  
  if(!missing(covariate)) {
    if(class(covariate) == "formula") {
      if(!missing(subset)) dat <- mcpl.prepare(formula=covariate,data=data,subset=subset,base=base) else
      dat <- mcpl.prepare(formula=covariate,data=data,base=base)
      z <- dat$x
    }
    else {
      z <- covariate
      if(is.matrix(covariate)) {
        if(NROW(z) != nrow(y)) stop("number of rows is not identical for responses and covariates")
      } else {
        if(length(z) != nrow(y)) stop("number of rows is not identical for responses and covariates") else {
          z <- as.matrix(z)
        }
      }
    }
  }
  
  # data <- as.data.frame(mf[,colnames(mf)!=learnpredname])
  
  parfill <- function(parameters) {
    if(family$family=="gaussian") if(is.null(parameters$sd)) parameters$sd <- 1 else parameters$sd <- parameters$sd
    if(is.null(parameters$coefficients)) parameters$coefficients <- rep(0,NCOL(x))
    parameters
  }
  
  if(is.null(ntimes) | replicate) {
    # intialize
    parameters <- parfill(parameters)
  } else {
    # setup a parlist
    nrep <- length(ntimes)
    # check structure of supplied list
    if(all(lapply(parameters,is.list)) && length(parameters)==nrep) {
      for(i in 1:nrep) {
        parameters[[i]] <- parfill(parameters[[i]])
      }
    } else {
      parameters <- parfill(parameters)
      parameters <- rep(list(parameters),nrep)
    }
  }

  if(missing(parStruct)) {
    tfix <- NULL
    if(!missing(fixed)) tfix <- fixed
    parStruct <- ParStruct(parameters,replicate=replicate,
                    fixed=tfix,ntimes=ntimes)
  }

  if(is.null(ntimes)) ntimes <- nrow(y)
  nTimes <- nTimes(ntimes)

#  if(is.null(learnpredname)) {
#    covariate <- FALSE
#    learnpredname <- ""
#  } else {
#    covariate <- TRUE
#  }
  mod <- new("GlmResponse",
    x = x,
    y = y,
    z = if(exists("z")) z else matrix(ncol=0,nrow=0),
    parStruct=parStruct,
    nTimes=nTimes,
    family=family,
    #formula=formula,
    combine=combine
  )
  mod <- runm(mod)
  mod
}

Try the mcplR package in your browser

Any scripts or data that you put into this service are public.

mcplR documentation built on May 2, 2019, 4:42 p.m.