This example is caused by R's messing with formula environments.
That's explained in ?formula, but should it not be explained in
?model.frame ?
Simple test:
data <- data.frame(y=rnorm(100), x=1:100)
testit <- function(formula)
{
weights <- runif(100)
glm(formula, weights=weights, data=data)
}
testit(y ~ x)
weights is looked for in the environment of the formula, not of the call
to glm.
--
Brian D. Ripley, ripley@stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272860 (secr)
Oxford OX1 3TG, UK Fax: +44 1865 272595
---------- Forwarded message ----------
Date: Tue, 21 Aug 2001 15:47:34 +0100 (BST)
From: Prof Brian Ripley <ripley@stats.ox.ac.uk>
To: Isabelle ZABALZA <isabelle.zabalza-mezghani@ifp.fr>
Subject: Re: [R] Problem using GLM in a loop
It's a scope problem: use
data$weights <- 1/sig
and it will work. It's finding the system function weights, not your
variable.
I don't really understand this one as yet.
On Tue, 21 Aug 2001, Isabelle ZABALZA wrote:
Hello,
I am try to perform a modeling which is relevant in a strongly
heteroscedastic context.
So I perform a dual modeling (modeling of both mean and variance of a
response) in using the following loop:
jointmod <- function(formula, data, itercrit=10,devcrit=0.0001)
{
#
# Init step
#
init <- glm(formula=formula,family=gaussian, data=data)
response <- resid(init,type="response") + predict(init,type="response")
mu <- predict(init,type="response")
iter <- 1
dev <- init$deviance
if (dev != 0) vardev <- 1
else stop(message="Null deviance in initialisation step")
endformula_ strsplit(formula,"~")[[3]]
formulavar_paste("d","~")
formulavar_as.formula(paste(formulavar,endformula))
#
# Modeling loop
#
while ((iter <= itercrit)| (vardev > devcrit)) {
d_(response - mu)^2
modvar_glm(formula=formulavar,family=Gamma(link=log),data=data)
sig <- predict(modvar,type="response")
weights <- 1/sig
modesp <- glm(formula=formula,
family=gaussian,data=data,weights=weights)
mu <- predict(modesp,type="response")
devold <- dev
dev <- modesp$deviance
vardev <- (devold -dev)/devold
iter <- iter + 1
}
list(modesp=modesp, modvar=modvar, iter=iter, vardev=vardev)
}
This program evaluation always stops when evaluating the "modesp" line :
modesp <- glm(formula=formula,
family=gaussian,data=data,weights=weights)
with the following error:
error in
model.frame(formula,rownames,variables,varnames,extras,extranames, :
invalid variable type.
The same program without specifying weights runs.
I don't understand this behavior since when I execute the glm command
with weights out of the program, it runs.
I run R1.3.0 under Windows.
If someone can help me ...
Isabelle Zabalza-Mezghani
--
Isabelle Zabalza-Mezghani Tel : 01 47 52 61 99
Institut Français du Pétrole E-mail : isabelle.zabalza-mezghani@ifp.fr
1-4 Av. Bois Preau - Bat Lauriers
92852 Rueil Malmaison Cedex, France
Brian D. Ripley, ripley@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._