root/MGET/Trunk/PythonPackage/src/GeoEco/Statistics/FitLMEForDataframe.r @ 981

# FitLMEForDataframe.r
#
# Copyright (C) 2011 Jason J. Roberts
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License (available in the file LICENSE.TXT)
# for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

FitLMEForDataframe <- function(fixed, data, outputModelFile, random, correlation="NULL", method=NULL, xVar=NULL, yVar=NULL, zVar=NULL, mVar=NULL, coordinateSystem=NULL, writeSummaryFile=TRUE, writeDiagnosticPlots=TRUE, plotFileFormat="png", res=1000.0, width=3000.0, height=3000.0, pointSize=10.0, bg="white")
{
    library(stats)
    
    # Fit the model.

    rPackage <- "nlme"
    library(rPackage, character.only=TRUE)
    message("Fitting the linear mixed model...")
    model <- lme(fixed, data=na.omit(data), random=random, correlation=correlation, method=method, keep.data=TRUE)
        
    # Print a summary of the model.
    
    modelSummary <- paste("\n",
                          "MODEL SUMMARY:\n",
                          "==============\n",
                          SummarizeModel(model), "\n",
                          sep="")

    writeLines(modelSummary)
    message("")

    # Write the output file.
    
    save(model, rPackage, fixed, random, xVar, yVar, zVar, mVar, coordinateSystem, file=outputModelFile, compress=FALSE)

    # Strip the extension, if one exists, off the output file path.
    
    outputFilePrefix <- file.path(dirname(outputModelFile), strsplit(basename(outputModelFile), ".", fixed=TRUE)[[1]])[1]

    # Write the summary file, if requested.
    
    if (writeSummaryFile)
    {
        f <- file(sprintf("%s_summary.txt", outputFilePrefix), "wt")
        tryCatch(writeLines(modelSummary, f), finally=close(f))
    }
    
    # Do all of the plotting in a tryCatch block, so that we can turn
    # off open graphics devices in the event of an error.

    tryCatch(
    {
        if (writeDiagnosticPlots)
        {
            # Plot the standardized residuals vs. fitted values.
            
##            if (plotFileFormat == "emf")
##                win.metafile(sprintf("%s_residuals_vs_fitted.emf", outputFilePrefix), width=width, height=height, pointsize=pointSize)
##            else
##                png(sprintf("%s_residuals_vs_fitted.png", outputFilePrefix), res=res, width=width, height=height, pointsize=pointSize, bg=bg)

##            model <<- model
##            plot.lme(model, resid(., type = "p") ~ fitted(.))
            
##            dev.off()
        }
    }, finally = {
        if (dev.cur() != 1)
            dev.off()
    })
}