SCM

SCM Repository

[matrix] Diff of /pkg/R/lmer.R
ViewVC logotype

Diff of /pkg/R/lmer.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 810, Thu Jul 14 13:48:38 2005 UTC revision 811, Thu Jul 14 15:37:41 2005 UTC
# Line 810  Line 810 
810            print(reMat, quote = FALSE)            print(reMat, quote = FALSE)
811        })        })
812    
813  glmmMCMC <- function(obj, nsamp = 1, alpha = 1, beta = 1, burnIn =  glmmMCMC <- function(obj, nsamp = 1, alpha = 1, beta = 1, saveb = FALSE, verbose = FALSE)
                      100, thinning = 5, b = FALSE, verbose = FALSE)  
814  {  {
     ## Perform N cycles of the chain.  This function is defined here  
     ## so it has access to objects in the environment  
     Ncycles <- function(state, N) {  
         fixed <- state$fixed  
         varc <- state$varc  
         b <- state$b  
         for (i in seq(len = N)) {  
             ## sample from the conditional distribution of beta given b and y  
             fixed <- .Call("glmer_fixed_update", GSpt, b,  
                            fixed, PACKAGE = "Matrix")  
             ## sample from the conditional distribution of b given beta, varc and y.  
             b <- .Call("glmer_ranef_update", GSpt, fixed, varc,  
                        b, PACKAGE = "Matrix")  
             ## sample from the conditional distribution of varc given b  
             varc <- 1/rgamma(1, shape = shape,  
                              scale = 1/(sum(b[[1]]^2)/2 + betainv))  
         }  
         list(fixed = fixed, varc = varc, b = b)  
     }  
   
815      ## Check arguments      ## Check arguments
816      if (!inherits(obj, "lmer")) stop("obj must be of class `lmer'")      if (!inherits(obj, "lmer")) stop("obj must be of class `lmer'")
817      if (obj@family$family == "gaussian" && obj@family$link == "identity")      if (obj@family$family == "gaussian" && obj@family$link == "identity")
818          warn("glmmMCMC not indended for Gaussian family with identity link")          warn("glmmMCMC not indended for Gaussian family with identity link")
819      if (length(obj@Omega) > 1 || obj@nc[1] > 1)      if (length(obj@Omega) > 1 || obj@nc[1] > 1)
820          stop("glmmMCMC currently defined for models with a single variance component")          stop("glmmMCMC currently defined for models with a single variance component")
821      cv <- lmerControl()      cv <- Matrix:::lmerControl()
822      if (verbose) cv$msVerbose <- 1      if (verbose) cv$msVerbose <- 1
823      family <- obj@family      family <- obj@family
824      frm <- obj@frame      frm <- obj@frame
# Line 867  Line 846 
846                 .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))                 .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
847      mer <- as(obj, "mer")      mer <- as(obj, "mer")
848    
849      ## establish the GS object      ## establish the GS object and the ans matrix
850      eta <- glm.fit$linear.predictors # perhaps later change this to obj@fitted?      eta <- glm.fit$linear.predictors # perhaps later change this to obj@fitted?
851      wts <- glm.fit$prior.weights      wts <- glm.fit$prior.weights
852      wtssqr <- wts * wts      wtssqr <- wts * wts
# Line 885  Line 864 
864      ## FIXME: The definition of shape is not general      ## FIXME: The definition of shape is not general
865      shape <- length(levels(mer@flist[[1]]))/2 + alpha      shape <- length(levels(mer@flist[[1]]))/2 + alpha
866      betainv <- 1/beta      betainv <- 1/beta
867        fixed <- obj@fixed
868        varc <- .Call("lmer_coef", mer, 2, PACKAGE = "Matrix")
869        b <- .Call("lmer_ranef", mer, PACKAGE = "Matrix")
870        ans <- array(0, c(nsamp, length(fixed) + length(varc)),
871                     list(NULL,
872                          c(names(fixed),
873                            paste("varc", seq(along = varc), sep = ''))))
874        if (saveb) {
875            blen <- length(unlist(b, recursive = TRUE))
876            ans <- cbind(ans, array(0, c(nsamp, blen),
877                                    list(NULL,
878                                         paste("b", 1:blen, sep = ''))))
879        }
880    
881      ## Perform burnIn cycles and establish the results matrix      ## create the samples
     state <- Ncycles(list(fixed = obj@fixed,  
                           varc = .Call("lmer_coef", mer, 2,  
                           PACKAGE = "Matrix"),  
                           b = .Call("lmer_ranef", mer,  
                           PACKAGE = "Matrix")),  
                      burnIn)  
     ans <- matrix(0, nrow = nsamp,  
                   ncol = length(state$fixed) +  
                   length(state$varc) +  
                   if (b) length(unlist(state$b,  
                                        recursive = TRUE)) else 0)  
882      for (i in 1:nsamp) {      for (i in 1:nsamp) {
883          state <- Ncycles(state, thinning)          ## sample from the conditional distribution of beta given b and y
884          ans[i,] <-          fixed <- .Call("glmer_fixed_update", GSpt, b,
885              c(state$fixed, state$varc,                             fixed, PACKAGE = "Matrix")
886                if (b) unlist(state$b, recursive = TRUE) else numeric(0))          ## sample from the conditional distribution of b given beta, varc and y.
887            b <- .Call("glmer_ranef_update", GSpt, fixed, varc,
888                       b, PACKAGE = "Matrix")
889            ## sample from the conditional distribution of varc given b
890            varc <- 1/rgamma(1, shape = shape,
891                             scale = 1/(sum(b[[1]]^2)/2 + betainv))
892            if (saveb)
893                ans[i,] <- c(fixed, varc, unlist(b, recursive = TRUE))
894            else
895                ans[i,] <- c(fixed, varc)
896      }      }
897      class(ans) <- "mcmc"      class(ans) <- "mcmc"
898      ans      ans

Legend:
Removed from v.810  
changed lines
  Added in v.811

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge