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 at ifp.fr
1-4 Av. Bois Preau - Bat Lauriers
92852 Rueil Malmaison Cedex, France
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
https://stat.ethz.ch/pipermail/r-help/attachments/20010821/ead5b473/attachment.html
I have done something similar, following the method described in Aitkin (1997, Applied Statistics, 36:332-339), and it worked fine. I am sending you my code (privately, it is very untidy :-( ) hoping this will help. Emmanuel Paradis At 16:23 21/08/01 +0200, Isabelle ZABALZA <isabelle.zabalza-mezghani at ifp.fr> 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 Franais du Ptrole E-mail : isabelle.zabalza-mezghani at ifp.fr 1-4 Av. Bois Preau - Bat Lauriers 92852 Rueil Malmaison Cedex, France -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help 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-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Isabelle ZABALZA <isabelle.zabalza-mezghani at ifp.fr> writes:> 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.This has to do with a rather awkward environment issue, which I'm not sure we have solved correctly. Your problem is really not the loop, but the fact that you are using glm inside a function which is passed the formula from another environment. Variables in a formula are taken from the supplied data frame and if not found in the environment of the formula. This is so that it would work to say e.g. y^lambda ~ x, where lambda is a free variable. However, the same treatment applies to special arguments like weight and offset, i.e. they are treated as names and not subjected to the usual parameter passing conventions. So the weights variable in the current environment is never found. Probably the easiest way out is to set environment(formula) <- environment() which might kill you if you're actually using a free variable as indicated, but chances are that you're not. Alternatively, but more than a bit icky, assign "weights" in the formula environment. -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help 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-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._