Rutils/maybe-not-useful/rsqu.r

In manfredo89/ED2io: Manages Input/Output for Ecosystem Demography 2 Model

#==========================================================================================#
#==========================================================================================#
#      This function computes the R2 for linear model objects.  Currently this has been    #
# tested with lm, glm, and rlm.                                                            #
#------------------------------------------------------------------------------------------#
rsqu <<- function(object){

   #------ Get original dependent variable. -----------------------------------------------#
   ymodel     = object$model[,1]
   ymodel.bar = mean(ymodel)
   #---------------------------------------------------------------------------------------#


   #---------------------------------------------------------------------------------------#
   #      Get the number of data points and degrees of freedom.                            #
   #---------------------------------------------------------------------------------------#
   ny     = length(ymodel)
   np     = length(object$coefficients)
   df.tot = ny - 1
   df.exp = ny - np - 1
   #---------------------------------------------------------------------------------------#


   #----- Get the total sum of squares. ---------------------------------------------------#
   ssq.tot  = sum2(ymodel-ymodel.bar)^2
   var.tot  = ssq.tot / df.tot
   #---------------------------------------------------------------------------------------#


   #----- Get the residual sum of squares. ------------------------------------------------#
   ssq.res  = sum2(object$residuals)^2
   var.res  = ssq.res / df.exp
   #---------------------------------------------------------------------------------------#


   #---------------------------------------------------------------------------------------#
   #      Result contains all statistics used to define R2 and adjusted R2.                #
   #---------------------------------------------------------------------------------------#
   ans = c( n         = ny
          , p         = np
          , df.tot    = df.tot
          , df.exp    = df.exp
          , sd.tot    = sqrt(var.tot)
          , sd.res    = sqrt(var.res)
          , r.squared = 1.0 - ssq.res / ssq.tot
          , adj.r2    = 1.0 - var.res / var.tot
          )#end c
   #---------------------------------------------------------------------------------------#

   return(ans)
}#end rsqu
#==========================================================================================#
#==========================================================================================#




#==========================================================================================#
#==========================================================================================#
#     Similar to rsqu, but applicable to other types of models.                            #
#------------------------------------------------------------------------------------------#
rsqu.gen <<- function(yobs,yhat,np=1,na.rm=TRUE,adjusted=TRUE){

   #------ Discard data in case na.rm = TRUE. ---------------------------------------------#
   if (na.rm){
      keep = is.finite(yobs) & is.finite(yhat)
      yobs = yobs[keep]
      yhat = yhat[keep]
   }#end if (na.rm)
   yres = yobs - yhat
   #---------------------------------------------------------------------------------------#



   #---------------------------------------------------------------------------------------#
   #     In case length is less than or equal to np, return NA.                            #
   #---------------------------------------------------------------------------------------#
   ny = length(yobs)
   if (ny <= (np+adjusted)){
      ans = NA
      return(ans)
   }#end if
   #---------------------------------------------------------------------------------------#



   #------ Get original dependent variable. -----------------------------------------------#
   yobs.bar = mean(yobs)
   #---------------------------------------------------------------------------------------#


   #---------------------------------------------------------------------------------------#
   #      Get the number of data points and degrees of freedom.                            #
   #---------------------------------------------------------------------------------------#
   df.tot = ny - 1
   df.exp = ny - np
   #---------------------------------------------------------------------------------------#


   #----- Get the total sum of squares. ---------------------------------------------------#
   ssq.tot  = sum((yobs-yobs.bar)^2)
   var.tot  = ssq.tot / df.tot
   #---------------------------------------------------------------------------------------#


   #----- Get the residual sum of squares. ------------------------------------------------#
   ssq.res  = sum(yres^2)
   var.res  = ssq.res / df.exp
   #---------------------------------------------------------------------------------------#



   #---------------------------------------------------------------------------------------#
   #      Result contains all statistics used to define R2 and adjusted R2.                #
   #---------------------------------------------------------------------------------------#
   if (adjusted){
      ans = 1.0 - var.res / var.tot
   }else{
      ans = 1.0 - ssq.res / ssq.tot
   }#end if (adjusted)
   #---------------------------------------------------------------------------------------#

   return(ans)
}#end rsqu.gen
#==========================================================================================#
#==========================================================================================#