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 706, Fri Apr 29 15:42:12 2005 UTC revision 707, Mon May 2 17:27:53 2005 UTC
# Line 282  Line 282 
282                invisible(object)                invisible(object)
283            })            })
284    
285    
286  setMethod("lmer", signature(formula = "formula"),  setMethod("lmer", signature(formula = "formula"),
287            function(formula, family, data,            function(formula, family, data,
288                     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),                     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
# Line 299  Line 300 
300                                  'for a generalized linear mixed model.',                                  'for a generalized linear mixed model.',
301                                  '\nUsing method = "PQL".\n'))                                  '\nUsing method = "PQL".\n'))
302            }            }
303            if (method %in% c("Laplace", "AGQ"))  ##           if (method %in% c("Laplace", "AGQ"))
304                stop("Laplace and AGQ methods not yet implemented")  ##               stop("Laplace and AGQ methods not yet implemented")
305              if (method %in% c("AGQ"))
306                  stop("AGQ method not yet implemented")
307            gVerb <- getOption("verbose")            gVerb <- getOption("verbose")
308                                          # match and check parameters                                          # match and check parameters
309            controlvals <- do.call("lmerControl", control)            controlvals <- do.call("lmerControl", control)
# Line 418  Line 421 
421            controlvals$msMaxIter <- msMaxIter.orig            controlvals$msMaxIter <- msMaxIter.orig
422    
423    
424            if (FALSE) ## Laplace  ###          if (TRUE) ## Laplace
425            {  ###          {
426                ## Need to optimize L(theta, beta) using Laplace approximation                ## Need to optimize L(theta, beta) using Laplace approximation
427    
428                ## Things needed for that:                ## Things needed for that:
# Line 553  Line 556 
556                        Omega[[i]] <- Omega[[i]] + t(Omega[[i]])                        Omega[[i]] <- Omega[[i]] + t(Omega[[i]])
557                        diag(Omega[[i]]) <- diag(Omega[[i]]) / 2                        diag(Omega[[i]]) <- diag(Omega[[i]]) / 2
558    
559                        ## want log of `const det(Omega) exp(-1/2 b' Omega b )`                    ## want log of `const det(Omega) exp(-1/2 b'
560                        ## i.e., const + log det(Omega) - .5 * (b' Omega b)                    ## Omega b )` i.e., const + log det(Omega) - .5
561                        ## FIXME: need to adjust for sigma^2 for appropriate models (easy)                    ## * (b' Omega b)
562                        ## these are all the b'Omega b, summed as they eventually need to be  
563                        ## think of this as sum(rowSums((ranefs[[i]] %*% Omega[[i]]) * ranefs[[i]]))                    ## FIXME: need to adjust for sigma^2 for
564                      ## appropriate models (easy).  These are all the
565                      ## b'Omega b, summed as they eventually need to
566                      ## be.  Think of this as
567                      ## sum(rowSums((ranefs[[i]] %*% Omega[[i]]) *
568                      ## ranefs[[i]]))
569    
570                        ranef.loglik.det <- nrow(ranefs[[i]]) *                        ranef.loglik.det <- nrow(ranefs[[i]]) *
571                            determinant(Omega[[i]], logarithm = TRUE)$modulus/2                            determinant(Omega[[i]], logarithm = TRUE)$modulus/2
# Line 570  Line 578 
578                        log.jacobian <-                        log.jacobian <-
579                            sum(log(abs(apply(reducedObj@bVar[[i]],                            sum(log(abs(apply(reducedObj@bVar[[i]],
580                                              3,                                              3,
581                                              function(x) sum(diag(x)))  
582                                            ## next line depends on
583                                            ## whether bVars are
584                                            ## variances or Cholesly
585                                            ## factors
586    
587                                            ## function(x) sum(diag(x)))
588                                            function(x) sum(diag( La.chol( x ) )))
589                                        )))                                        )))
590    
591                        ## the constant terms from the r.e. and the final                        ## the constant terms from the r.e. and the final
# Line 594  Line 609 
609    
610                if (method == "Laplace")                if (method == "Laplace")
611                {                {
612                  ##Rprof() # trying to figure out if C-ifying bhat is worthwhile
613                    ## no analytic gradients or hessians                    ## no analytic gradients or hessians
614                    optimRes <-                    optimRes <-
615                        optim(fn = devLaplace,                        optim(fn = devLaplace,
616                              par = c(fixef(obj), coef(obj, unconst = TRUE)),                          par =
617                            c(fixef(obj),
618                              .Call("lmer_coef",
619                                    obj,
620                                    TRUE,
621                                    PACKAGE = "Matrix")),
622                            ## WAS: coef(obj, unconst = TRUE)),
623                              method = "BFGS", hessian = TRUE,                              method = "BFGS", hessian = TRUE,
624                              control = list(trace = getOption("verbose"),                              control = list(trace = getOption("verbose"),
625                              reltol = controlvals$msTol,                              reltol = controlvals$msTol,
# Line 615  Line 637 
637                        print(fixef(obj))                        print(fixef(obj))
638                        print(optimRes$par[seq(length = responseIndex - 1)])                        print(optimRes$par[seq(length = responseIndex - 1)])
639                        cat("(Unconstrained) variance coefficients:\n")                        cat("(Unconstrained) variance coefficients:\n")
640                        print(coef(obj, unconst = TRUE))  
641                        coef(obj, unconst = TRUE) <-                    print(
642                            optimRes$par[responseIndex:length(optimRes$par)]                          .Call("lmer_coef",
643                        print(coef(obj, unconst = TRUE))                                obj,
644                                  TRUE,
645                                  PACKAGE = "Matrix"))
646    
647                      ##coef(obj, unconst = TRUE) <-
648                      ##    optimRes$par[responseIndex:length(optimRes$par)]
649                      ##print(coef(obj, unconst = TRUE))
650                      print( optimRes$par[responseIndex:length(optimRes$par)] )
651                    }                    }
652    
653                    ## need to calculate likelihood also need to store new                ## need to calculate likelihood.  also need to store
654                    ## estimates of fixed effects somewhere (probably cannot                ## new estimates of fixed effects somewhere
655                    ## update standard errors)                ## (probably cannot update standard errors)
656                                            #Rprof(NULL)
657                }                }
658                else                else
659                {                {
# Line 644  Line 674 
674                ## have the 'correct' random effects in reducedObj.                ## have the 'correct' random effects in reducedObj.
675    
676                loglik <- devLaplace(optpars)                loglik <- devLaplace(optpars)
677                print(loglik)            ##print(loglik)
678                ff <- optpars[1:(responseIndex-1)]                ff <- optpars[1:(responseIndex-1)]
679                names(ff) <- names(fixef(obj))                names(ff) <- names(fixef(obj))
680    
681                if (!x) mmats <- list()                if (!x) mmats <- list()
682    
683            }  ###      }
684    
685            obj            obj
686        })        })
687    
688    
689  ## calculates degrees of freedom for fixed effects Wald tests  ## calculates degrees of freedom for fixed effects Wald tests
690  ## This is a placeholder.  The answers are generally wrong.  It will  ## This is a placeholder.  The answers are generally wrong.  It will
691  ## be very tricky to decide what a 'right' answer should be with  ## be very tricky to decide what a 'right' answer should be with

Legend:
Removed from v.706  
changed lines
  Added in v.707

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