[R] Problem using GLM in a loop

[Emmanuel Paradis]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._