[Rd] model.matrix evaluation challenges
Ben Bolker
bolker at ufl.edu
Tue Aug 11 16:22:31 CEST 2009
that seems to work, thank you.
would anyone care to explain *why* it works ... ? or where
I could go to read more about why it works ... ?
cheers
Ben Bolker
Felix Andrews wrote:
> how about...
>
> nrow(with(cc, model.matrix(params, data=environment())))
>
> cheers
> -Felix
>
> 2009/8/10 Ben Bolker <bolker at ufl.edu>:
>> I am having difficulty with evaluation/environment construction
>> for a formula to be evaluated by model.matrix(). Basically, I
>> want to construct a model matrix that first looks in "newdata"
>> for the values of the model parameters, then in "object at data".
>> Here's what I've tried:
>>
>> 1. model.matrix(~f,data=c(newdata,object at data)) -- fails because
>> something (terms()?) within model.matrix() tries to coerce "data"
>> to a data frame -- but newdata and object at data have different
>> numbers of rows.
>>
>> 2. with(c(newdata,object at data),model.matrix(~f)) works --
>> BUT not when ~f is assigned to another object (as it will be
>> when it is being dug out of a previously fit model).
>>
>> 3. If "f" exists as a column within newdata and/or object at data,
>> but the formula is assigned to another object (see below), then
>> with(c(newdata,object at data),model.matrix(z)) fails -- because
>> model.matrix() is explicitly evaluating the variables in the formula
>> in the environment of z (i.e., ignoring the first argument of "with" ...)
>>
>>
>> Any advice on how to solve this without making a bigger mess?
>>
>> sincerely
>> Ben Bolker
>>
>>
>> ## set up a data frame for prediction
>>
>> set.seed(1001)
>> f = factor(rep(letters[1:4],each=20))
>> x = runif(80)
>> u = rnorm(4)
>> y = rnorm(80,mean=2+x*(3+u[f]),sd=0.1)
>> dat = data.frame(f,x,y)
>>
>> ## fit a model ... could easily do by lm() but want to
>> ## demonstrate the problem
>>
>> library(bbmle)
>> m1 = mle2(y~dnorm(a+b*x,sd=exp(logs)),parameters=list(b~f),data=dat,
>> start=list(a=0,b=2,logs=-3))
>>
>> ## data frame for prediction
>> pp0 = expand.grid(x=seq(0,1,length=11),
>> f=levels(dat$f))
>>
>> ## combine frame and model data: have to keep the model data
>> ## around, because it contain other information needed for
>> ## prediction.
>>
>> cc = c(pp0,m1 at data)
>> try(mm <- model.matrix(~f,cc)) ## fails -- different number of rows
>> ## (at some point R tries to coerce cc to a data frame)
>> nrow(with(cc,model.matrix(~f))) ## works ... but ...
>>
>> ## here's what happens within predict.mle2()
>> ## extract information about parameter models
>> params <- eval(m1 at call$parameters)
>> params <- params[[1]]
>> params[[2]] <- NULL
>>
>> nrow(with(cc,model.matrix(params)))
>> rm(f)
>> ## object "f" not found -- because model.matrix
>> ## evaluates in the parent environment of params ...
>> try(nrow(with(cc,model.matrix(params))))
>>
>>
>>
>> --
>> Ben Bolker
>> Associate professor, Biology Dep't, Univ. of Florida
>> bolker at ufl.edu / www.zoology.ufl.edu/bolker
>> GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc
>>
>>
>> ______________________________________________
>> R-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>>
>
>
>
--
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / www.zoology.ufl.edu/bolker
GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 261 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20090811/5705b1d1/attachment.bin>
More information about the R-devel
mailing list