Re: [Rd] Proposal: model.data

From: Paul Johnson <pauljohn32_at_gmail.com>
Date: Fri, 04 May 2012 11:25:58 -0500

On Thu, May 3, 2012 at 12:19 PM, Brian G. Peterson <brian_at_braverock.com> wrote:
> On Thu, 2012-05-03 at 12:09 -0500, Paul Johnson wrote:

>> Greetings:
>>
>> On Thu, May 3, 2012 at 11:36 AM, Brian G. Peterson <brian_at_braverock.com> wrote:
>> > On Thu, 2012-05-03 at 10:51 -0500, Paul Johnson wrote:
>> >> If somebody in R Core would like this and think about putting it, or
>> >> something like it, into the base, then many chores involving predicted
>> >> values would become much easier.
>> >>
>> > Why does this need to be in base?  Implement it in a package.
>> >>

So, nobody agrees with me that R base should have model.data? Too bad! You could save us a lot of effort. I've found different efforts to get the same work done in termplot and all of the implementations of predict. And just about every regression related package has its own approach.

If there were one good way that would always work, just think how convenient it would be for all those function & package writers. Oh, well. I'm not saying that my model.data is perfect, just that I wish there were a perfect one :)

Yesterday, I realized that predict.nls probably has to deal with this problem so I studied that and have yet another version of model.data to propose to you. I'm using this in my regression-support package rockchalk, so you don't need to give me Brian Peterson's advice to "put it in a package".

The idea here is to take variables from a fitted model's data if it can find them, and then grab variables from the parent environment IF they have the correct length. This means we ignore variables like d in poly(x, d) because the variable d is not of the same length as the variables in the model.

##' Creates a "raw" (UNTRANSFORMED) data frame equivalent
##' to the input data that would be required to fit the given model.
##'
##' Unlike model.frame and model.matrix, this does not return transformed
##' variables.
##'
##' @param model A fitted regression model in which the data argument
##' is specified. This function will fail if the model was not fit
##' with the data option.
##' @return A data frame
##' @export
##' @author Paul E. Johnson <pauljohn@@ku.edu>
##' @example inst/examples/model.data-ex.R
model.data <- function(model){

    #from nls, returns -1 for missing variables     lenVar <- function(var, data) tryCatch(length(eval(as.name(var),

                         data, env)), error = function(e) -1)
    fmla <- formula(model)
    varNames <- all.vars(fmla) ## all variable names
    ## varNames includes d in poly(x,d), possibly other "constants"
    ## varNamesRHS <- all.vars(formula(delete.response(terms(model))))
    ## previous same as nls way?

    fmla2 <- fmla
    fmla2[[2L]] <- 0
    varNamesRHS <- all.vars(fmla2)
    varNamesLHS <- setdiff(varNames, varNamesRHS)
    env <- environment(fmla)
    if (is.null(env))

        env <- parent.frame()

    dataOrig <- eval(model$call$data, environment(formula(model)))     rndataOrig <- row.names(dataOrig)
    n <- sapply(varNames, lenVar, data=dataOrig)     targetLength <- length(eval(as.name(varNamesLHS[1]), dataOrig, env))     varNames <- varNames[ n == targetLength ]     ldata <- lapply(varNames, function(x) eval(as.name(x), dataOrig, env))     names(ldata) <- varNames
    data <- data.frame(ldata[varNames])
    if (!is.null(rndataOrig)) row.names(data) <- rndataOrig     ## remove rows listed in model's na.action     ## TODO: question: what else besides OMIT might be called for?     if ( !is.null(model$na.action)){

        data <- data[ -as.vector(model$na.action), ]     }
    ## keep varNamesRHS that exist in datOrig     attr(data, "varNamesRHS") <- setdiff(colnames(data), varNamesLHS)     invisible(data)
}

And some example output:
> ## check if model.data works when there is no data argument
> set.seed(12345)
> x1 <- rpois(100, l=6)
> x2 <- rnorm(100, m=50, s=10)
> x3 <- rnorm(100)
> y <- rnorm(100)
> m0 <- lm(y ~ x1 + x2 + x3)
> m0.data <- model.data(m0)
> x1[4:5] <- NA
> m0 <- lm(y ~ x1 + x2 + x3)
> m0.data <- model.data(m0)
> head(m0.data)

           y x1 x2 x3

1 -0.8086741  7 44.59614 -1.6193283
2  1.0011198  9 69.47693  0.5483979
3  0.4560525  8 50.53590  0.1952822
6  0.6417692  4 52.77954 -0.2509466
7 -0.4150210  5 56.91171  1.6993467
8 -0.4595757  6 58.23795 -0.3442988

> x1 <- rpois(100, l=6)
> x2 <- rnorm(100, m=50, s=10)
> x3 <- rnorm(100)
> xcat1 <- gl(2,50, labels=c("M","F"))
> xcat2 <- cut(rnorm(100), breaks=c(-Inf, 0, 0.4, 0.9, 1, Inf), labels=c("R", "M", "D", "P", "G"))
> dat <- data.frame(x1, x2, x3, xcat1, xcat2)
> rm(x1, x2, x3, xcat1, xcat2)
> xcat1n <- with(dat, contrasts(xcat1)[xcat1, ,drop=FALSE])
> xcat2n <- with(dat, contrasts(xcat2)[xcat2, ])
> STDE <- 20
> dat$y <- 0.03 + 0.8*dat$x1 + 0.1*dat$x2 + 0.7*dat$x3 + xcat1n %*% c(2) + xcat2n %*% c(0.1,-2,0.3, 0.1) + STDE*rnorm(100)
> rownames(dat$y) <- NULL
> m1 <- lm(y ~ poly(x1, 2), data=dat)
> m1.data <- model.data(m1)
> head(m1.data)

           y x1

1  56.336279  7
2  47.823205  5
3   9.296108  6
4  16.213508  5
5 -16.922331  3
6  10.639724  7

> attr(m1.data, "varNamesRHS")

[1] "x1"
> d <- 2
> m2 <- lm(y ~ poly(x1, d), data=dat)
> m2.data <- model.data(m2)
> head(m2.data)

           y x1

1  56.336279  7
2  47.823205  5
3   9.296108  6
4  16.213508  5
5 -16.922331  3
6  10.639724  7

> attr(m2.data, "varNamesRHS")

[1] "x1"
> m3 <- lm(y ~ log(10 + x1) + poly(x1, d) + sin(x2), data=dat)
> m3.data <- model.data(m3)
> head(m3.data)

           y x1 x2

1  56.336279  7 56.27965
2  47.823205  5 50.02144
3   9.296108  6 52.84378
4  16.213508  5 39.98221
5 -16.922331  3 43.82778
6  10.639724  7 58.28194

> attr(m3.data, "varNamesRHS")

[1] "x1" "x2"
> m4 <- lm(log(50+y) ~ log(d+10+x1) + poly(x1, 2), data=dat)
> m4.data <- model.data(m4)
> head(m4.data)

           y x1

1  56.336279  7
2  47.823205  5
3   9.296108  6
4  16.213508  5
5 -16.922331  3
6  10.639724  7

> attr(m4.data, "varNamesRHS")

[1] "x1"
> m4 <- lm(y ~ x1*x1, data=dat)
> m4.data <- model.data(m4)
> head(m4.data)

           y x1

1  56.336279  7
2  47.823205  5
3   9.296108  6
4  16.213508  5
5 -16.922331  3
6  10.639724  7

> attr(m4.data, "varNamesRHS")

[1] "x1"
> m4 <- lm(y ~ x1 + I(x1^2), data=dat)
> m4.data <- model.data(m4)
> head(m4.data)

           y x1

1  56.336279  7
2  47.823205  5
3   9.296108  6
4  16.213508  5
5 -16.922331  3
6  10.639724  7

> attr(m4.data, "varNamesRHS")

[1] "x1"
> dat$x1[sample(100, 5)] <- NA
> dat$y[sample(100, 5)] <- NA
> dat$x2[sample(100, 5)] <- NA
> dat$x3[sample(100,10)] <- NA
> m1 <- lm(y ~ log(10 + x1), data=dat)
> m1.data <- model.data(m1)
> head(m1.data)

           y x1

1  56.336279  7
2  47.823205  5
3   9.296108  6
4  16.213508  5
5 -16.922331  3
6  10.639724  7

> summarize(m1.data)

$numerics

         x1 y

0%    1.000 -31.9800
25%   5.000   0.9636
50%   6.000  13.8400
75%   7.000  27.5300
100% 12.000  71.1100
mean  5.933  14.8400
sd    2.336  20.4300
var   5.456 417.3000
NA's  0.000   0.0000
N    90.000  90.0000

$factors
NULL
> attr(m1.data, "varNamesRHS")

[1] "x1"
> m2 <- lm(y ~ log(x1 + 10), data=dat)
> m1.data <- model.data(m1)
> head(m1.data)

           y x1

1  56.336279  7
2  47.823205  5
3   9.296108  6
4  16.213508  5
5 -16.922331  3
6  10.639724  7

> summarize(m1.data)

$numerics

         x1 y

0%    1.000 -31.9800
25%   5.000   0.9636
50%   6.000  13.8400
75%   7.000  27.5300
100% 12.000  71.1100
mean  5.933  14.8400
sd    2.336  20.4300
var   5.456 417.3000
NA's  0.000   0.0000
N    90.000  90.0000

$factors
NULL
> attr(m1.data, "varNamesRHS")

[1] "x1"
> d <- 2
> m3 <- lm(log(50+y) ~ log(d+10+x1) + x2 + sin(x3), data=dat)
> m3.data <- model.data(m3)
> head(m3.data)

           y x1 x2 x3

2  47.823205  5 50.02144 -2.4669386
3   9.296108  6 52.84378  0.4847158
4  16.213508  5 39.98221 -0.9379723
5 -16.922331  3 43.82778  3.3307333
7  26.084587  3 49.15181  0.2204558
8   4.392061  8 45.65280  0.8762108

> summarize(m3.data)

$numerics

         x1 x2 x3 y

0%    1.000 27.630 -2.55600 -31.9800
25%   4.000 42.370 -0.52690   0.7905
50%   6.000 49.150  0.04162  11.5800
75%   7.000 54.390  0.71310  26.8900
100% 12.000 72.390  3.33100  71.1100
mean  5.883 48.800  0.02186  13.2300
sd    2.362  8.847  1.08900  20.0400
var   5.578 78.270  1.18700 401.8000
NA's  0.000  0.000  0.00000   0.0000
N    77.000 77.000 77.00000  77.0000

$factors
NULL
> attr(m3.data, "varNamesRHS")

[1] "x1" "x2" "x3"
> m4 <- lm(y ~ x1^3 + log(x2), data=dat)
> m4.data <- model.data(m4)
> summarize(m4.data)
$numerics

         x1 x2 y

0%    1.000 27.630 -31.9800
25%   4.000 42.630   0.8032
50%   6.000 49.710  12.2500
75%   7.000 54.910  27.2300
100% 12.000 72.390  71.1100
mean  5.918 49.220  14.0800
sd    2.372  8.808  20.0000
var   5.624 77.590 399.9000
NA's  0.000  0.000   0.0000
N    85.000 85.000  85.0000

$factors
NULL
> attr(m4.data, "varNamesRHS")

[1] "x1" "x2"
> m5 <- lm(y ~ x1 + I(x1^2) + cos(x2), data=dat)
> m5.data <- model.data(m5)
> head(m5.data)

           y x1 x2

1  56.336279  7 56.27965
2  47.823205  5 50.02144
3   9.296108  6 52.84378
4  16.213508  5 39.98221
5 -16.922331  3 43.82778
7  26.084587  3 49.15181

> summarize(m5.data)

$numerics

         x1 x2 y

0%    1.000 27.630 -31.9800
25%   4.000 42.630   0.8032
50%   6.000 49.710  12.2500
75%   7.000 54.910  27.2300
100% 12.000 72.390  71.1100
mean  5.918 49.220  14.0800
sd    2.372  8.808  20.0000
var   5.624 77.590 399.9000
NA's  0.000  0.000   0.0000
N    85.000 85.000  85.0000

$factors
NULL
> attr(m5.data, "varNamesRHS")

[1] "x1" "x2"
> x10 <- rnorm(100)
> x11 <- rnorm(100)
> m6 <- lm(y ~ x1 + I(x1^2) + cos(x2) + log(10 + x10) + sin(x11) + x10*x11, data=dat)
> m6.data <- model.data(m6)
> head(m6.data)

           y x1 x2 x10 x11

1  56.336279  7 56.27965 -0.4562360 -0.27005578
2  47.823205  5 50.02144 -1.4031737 -1.11355194
3   9.296108  6 52.84378 -0.2816067 -0.08319405
4  16.213508  5 39.98221  0.3353308  0.17284191
5 -16.922331  3 43.82778 -0.5184262  0.61754378
7  26.084587  3 49.15181  0.9037821  0.02170007

> dim(m6.data)

[1] 85 5
> summarize(m5.data)

$numerics

         x1 x2 y

0%    1.000 27.630 -31.9800
25%   4.000 42.630   0.8032
50%   6.000 49.710  12.2500
75%   7.000 54.910  27.2300
100% 12.000 72.390  71.1100
mean  5.918 49.220  14.0800
sd    2.372  8.808  20.0000
var   5.624 77.590 399.9000
NA's  0.000  0.000   0.0000
N    85.000 85.000  85.0000

$factors
NULL
> attr(m6.data, "varNamesRHS")

[1] "x1" "x2" "x10" "x11"
>

-- 
Paul E. Johnson
Professor, Political Science    Assoc. Director
1541 Lilac Lane, Room 504     Center for Research Methods
University of Kansas               University of Kansas
http://pj.freefaculty.org            http://quant.ku.edu

______________________________________________
R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Fri 04 May 2012 - 16:51:59 GMT

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

All messages

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 04 May 2012 - 17:21:07 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive