SH: response

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: doBy
-Version: 4.6.25.9004
+Version: 4.6.25.9016
 Title: Groupwise Statistics, LSmeans, Linear Estimates, Utilities 
 Authors@R: c(
     person(given = "Ulrich", family = "Halekoh",

NAMESPACE

@@ -69,6 +69,8 @@ S3method(tidy,linest_class)
 S3method(vcov,esticon_class)
 export(LE_matrix)
 export(LSmeans)
+export(add_pred)
+export(add_resid)
 export(aggregate_linest_list)
 export(as_lhs_chr)
 export(as_lhs_frm)
@@ -119,6 +121,7 @@ export(recodeVar)
 export(recode_var)
 export(recover_pca_data)
 export(renameCol)
+export(response)
 export(response_plot)
 export(rle2)
 export(sampleBy)
@@ -194,6 +197,8 @@ importFrom(stats,logLik)
 importFrom(stats,median)
 importFrom(stats,model.frame)
 importFrom(stats,model.matrix)
+importFrom(stats,model.response)
+importFrom(stats,model.weights)
 importFrom(stats,pchisq)
 importFrom(stats,pnorm)
 importFrom(stats,predict)
@@ -204,6 +209,7 @@ importFrom(stats,qt)
 importFrom(stats,resid)
 importFrom(stats,residuals)
 importFrom(stats,rstandard)
+importFrom(stats,rstudent)
 importFrom(stats,sd)
 importFrom(stats,step)
 importFrom(stats,summary.lm)

R/linear_modelling.R

@@ -1,11 +1,21 @@
 
 ##' @title Add predicted values of different types to dataframe
+##' 
 ##' @param data dataframe or tibble
 ##' @param model model object
 ##' @param var name of new variable in dataframe / tibble
 ##' @param type type of predicted value 
 ##' @return dataframe / tibble
 ##' @author Søren Højsgaard
+##' @examples
+##' data(cars)
+##' lm1 <- lm(dist ~ speed + I(speed^2), data=cars)
+##' lm1 |> response() |> head()
+##' cars <- cars |> add_pred(lm1)
+##' cars |> head()
+##' cars <- cars |> add_resid(lm1)
+##' cars
+##' 
 ##' @export
 add_pred <- function (data, model, var = "pred", type = NULL) 
 {
@@ -24,15 +34,25 @@ add_pred <- function (data, model, var = "pred", type = NULL)
 }
 
 ##' @title Add residuals of different types to dataframe
+##' 
 ##' @param data dataframe or tibble
 ##' @param model model object
 ##' @param var name of new variable in dataframe / tibble
 ##' @param type type of residual value 
 ##' @return dataframe / tibble
 ##' @author Søren Højsgaard
+##' @examples
+##' data(cars)
+##' lm1 <- lm(dist ~ speed + I(speed^2), data=cars)
+##' lm1 |> response() |> head()
+##' cars <- cars |> add_pred(lm1)
+##' cars |> head()
+##' cars <- cars |> add_resid(lm1)
+##' cars 
+##'
 ##' @export
-add_resid <- function (data, model, var = "resid", type) 
-{
+add_resid <- function (data, model, var = "resid", type) {
+
     resid2 <- function(model, type){
         UseMethod("resid2")
     }
@@ -49,15 +69,31 @@ add_resid <- function (data, model, var = "resid", type)
             return(stats::rstudent(model))                       
         }
 
-        return(stats::residuls(model))
+        return(stats::residuals(model))
     }
+    if (missing(type))
+        type="working"
+
+
     data[[var]] <- resid2(model, type)
+    data
 }
 
-##' @title Returns the response in a regression model.
-##' @param object Currently an lm og glm.
-##' @return A vector (or a matrix)
-##' @author Søren Højsgaard
+
+##' @title Get response variable from model
+##' @param object lm or glm object 
+##' @param data dataframe or tibble
+##' @param model model object
+##' @param var name of new variable in dataframe / tibble
+##' @param type type of residual value
+##' @examples
+##' data(cars)
+##' lm1 <- lm(dist ~ speed + I(speed^2), data=cars)
+##' lm1 |> response() |> head()
+##' cars <- cars |> add_pred(lm1)
+##' cars |> head()
+##' cars <- cars |> add_resid(lm1)
+##' cars 
 ##' @export
 response <- function(object){
 
@@ -75,7 +111,7 @@ response <- function(object){
     }
 
     obs.lm <- function(object) {
-        obs <- unname(model.response(model.frame(object)))        
+        obs <- model.response(model.frame(object))        
 
         if (is_glm(object)){
             wgt <- unname(model.weights(model.frame(object)))

R/plot_functions.r

@@ -4,37 +4,47 @@
 #' @param format The format of the plot (or a list of plots if format is "list")
 #'
 #' @examples
-#'
-#' m1 <- lm(speed ~ dist, data=cars)
+#' data(income)
+#' m1 <- lm(inc ~ race + educ, data=income)
 #' plot_lm(m1)
 #' plot_lm(m1, "2x2")
 #' plot_lm(m1, "1x4")
 #' plot_lm(m1, "4x1")
 #' plot_lm(m1, "list")
 #' 
 #' @export
-plot_lm <- function(lm_fit, format="2x2") {
-
+plot_lm <- function(lm_fit, format="2x2", global_aes=NULL) {
+    
     if(!inherits(lm_fit, "lm"))
         stop("'lm_fit' must inherit form 'lm'\n")
 
     format <- match.arg(format, c("2x2", "1x4", "4x1", "list"))
 
-    dd <- data.frame(fitted_values          = predict(lm_fit),
-                   ## residuals           = resid(lm_fit),
-                   observed_values        = predict(lm_fit)+resid(lm_fit),
-                   standardized_residuals = rstandard(lm_fit))
-
-    pl1 <- ggplot(dd, aes(x = .data$fitted_values, y = .data$standardized_residuals)) +
+    dd <- data.frame(
+        fitted_values          = predict(lm_fit),
+        resp        = predict(lm_fit)+resid(lm_fit),
+        std_res = rstandard(lm_fit))
+
+
+    pl1 <- ggplot(dd,
+                  aes(x = .data$fitted_values,
+                      y = .data$std_res)) +
         geom_point() + geom_hline(yintercept=c(0, -1.96, 1.95))
 
-    pl2 <- ggplot(dd, aes(x=.data$fitted_values, y=.data$observed_values)) +
+    pl2 <- ggplot(dd,
+                  aes(x=.data$fitted_values,
+                      y=.data$resp)) +
         geom_point() +  geom_smooth(method="lm", formula=y ~ x, se=FALSE)
 
-    pl3 <- ggplot(dd, aes(x=.data$fitted_values, y=sqrt(abs(.data$standardized_residuals)))) + geom_point()
+    pl3 <- ggplot(dd,
+                  aes(x=.data$fitted_values,
+                      y=sqrt(abs(.data$std_res)))) +
+        geom_point()
 
-    pl4 <- ggplot(dd, aes(sample = .data$standardized_residuals)) + stat_qq() + stat_qq_line() +
-        labs(y="standardized_residuals", x = "theoretical quantiles")
+    pl4 <- ggplot(dd,
+                  aes(sample = .data$std_res)) +
+        stat_qq() + stat_qq_line() +
+        labs(y="std_res", x = "theoretical quantiles")
 
     switch(format,
            "2x2"={ cowplot::plot_grid(pl1, pl2, pl3, pl4, nrow=2) },
@@ -46,6 +56,24 @@ plot_lm <- function(lm_fit, format="2x2") {
 }
 
 
+
+    ## aes_ <- c(list(x = predict(lm_fit),
+    ##                y = rstandard(lm_fit)),
+    ##           global_aes)
+
+    ## print(aes_)
+
+    ## aes_ <- ggpubr_create_aes(aes_)
+
+    ## print(aes_)
+
+    ## pl1 <- ggplot(dd, aes_) +
+        ## geom_point() + geom_hline(yintercept=c(0, -1.96, 1.95))
+
+    ## pl1 <- ggplot(dd,
+    ##               ggpubr_create_aes(list(x = .data$fitted_values,
+    ##                   y = .data$standardized_residuals))) +
+    ##     geom_point() + geom_hline(yintercept=c(0, -1.96, 1.95))