[R] Problem using GLM in a loop

[Isabelle ZABALZA]

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