[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