From pauljohn32 at gmail.com  Fri May  4 18:25:58 2012
From: pauljohn32 at gmail.com (Paul Johnson)
Date: Fri, 4 May 2012 11:25:58 -0500
Subject: [Rd] Proposal: model.data
In-Reply-To: <1336065561.29163.42.camel@brian-rcg>
References: <CAErODj8fUQW85O6hMzU40ACH3Nvzk7fBa2OnV=mx_MnNwqDC5A@mail.gmail.com>
	<1336063009.29163.39.camel@brian-rcg>
	<CAErODj8DhfYfn7y0kwWcUfKWCRVc-MZbEBTAuW_mEfQhnwfLSg@mail.gmail.com>
	<1336065561.29163.42.camel@brian-rcg>
Message-ID: <CAErODj8EH_Ht9LD7_UWtRsSHgGsaHhn3fZ2zOu4wnYsTk8eQ1Q@mail.gmail.com>

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



[Rd] Proposal: model.data
Paul Johnson