[R] Formula aruguments with NLS and model.frame()

[Spencer Graves]

      I haven't seen any replies to this post, so I will offer a couple 
of comments. 

      First, your example is entirely too complicated and not self 
contained for me to say much in the limited time I have for this.  I 
suggest you start by splitting your problem into several steps.  For 
example, will 'garchFit' allow 'formula.mean' to be something other than 
'~arma(p, q)', where p and q are nonnegative integers?  If no, I suggest 
you start by trying to produce your own modification to 'garchFit' so it 
will accept something like 'formula.mean=~z+arma(p, q)';  I suggest you 
give your local copy a different name like 'garchFitZ'.  Second, I 
suggest you cut your example data down to a minimum, e.g, 9 or 20 
observations, just enough so the algorithm won't die for some reason 
that would not occur with a larger data set but small enough so you can 
quickly print out and study every object the 'garchFit' algorithm 
produces. 

      Second, I suggest you use 'debug' to walk through your local 
version of 'garchFit' line by line.  I've found this to be a very 
powerful way to learn what happens in the internal environment of a 
function. 

      If you get stuck trying this, please submit another post including 
commented, minimal, self-contained, reproducible code, as suggested in 
the posting guide 'www.R-project.org/posting-guide.html'. 
      Hope this helps. 
      Spencer Graves

and provide commented, minimal, self-contained, reproducible code.



Joe W. Byers wrote:
> I could use some help understanding how nls parses the formula argument
> to a model.frame and estimates the model.  I am trying to utilize the
> functionality of the nls formula argument to modify garchFit() to handle
> other variables in the mean equation besides just an arma(u,v)
> specification.
>
> My nonlinear model is
>       y<-nls(t~a*sin(w*2*pi/365*id+p)+b*id+int,data=t1,
> 	start=list(w=.5,a=.1,p=.5,b=init.y$coef[2],int=init.y$coef[1] ),
> 	control=list(maxiter=1000000,minFactor=1e-18))
> where t is change in daily temperatures, id is just a time trend and the
> a*sin is a one year fourier series.
>
> I have tried to debug the nls code using the following code
> t1<-data.frame(t=as.vector(x),id=index(x))
> data=t1;
> formula <- as.formula(t ~ a *sin(w *2* pi/365 * id + p) + b * id + int);
>       varNames <- all.vars(formula)
>       algorithm<-'default';
>       mf <- match.call(definition=nls,expand.dots=FALSE,
>       call('nls',formula, data=parent.frame(),start,control = nls.control(),
>       algorithm = "default", trace = FALSE,
>       subset, weights, na.action, model = FALSE, lower = -Inf,
>       upper = Inf));
>       mWeights<-F;#missing(weights);
> 	start=list(w=.5,a=.1,p=.5,b=init.y$coef[2],int=init.y$coef[1] );
>       pnames <- names(start);
>        varNames <- varNames[is.na(match(varNames, pnames, nomatch = NA))]
>
> 	varIndex <- sapply(varNames,
> 		function(varName, data, respLength) {
>           	length(eval(as.name(varName), data))%%respLength == 0},
>           	 data, length(eval(formula[[2]], data))
>           );
> 	mf$formula <- as.formula(paste("~", paste(varNames[varIndex],
>           collapse = "+")), env = environment(formula));
> 	mf$start <- NULL;mf$control <- NULL;mf$algorithm <- NULL;
> 	mf$trace <- NULL;mf$model <- NULL;
>       mf$lower <- NULL;mf$upper <- NULL;
>       mf[[1]] <- as.name("model.frame");
>       mf<-evalq(mf,data);
>       n<-nrow(mf)
>       mf<-as.list(mf);
>       wts <- if (!mWeights)
>           model.weights(mf)
>       else rep(1, n)
>       if (any(wts < 0 | is.na(wts)))
>           stop("missing or negative weights not allowed")
>
>       m <- switch(algorithm,
>       		plinear = nlsModel.plinear(formula, mf, start, wts),
>       		port = nlsModel(formula, mf, start, wts, upper),
>       		nlsModel(formula, mf, start, wts));
>
> I am struggling with the environment issues associated with performing
> these operations.  I did not include the data because it is 9000 
> observations of temperature data.  If anyone would like the data, I can 
> provide it or a subset in a csv file.
>
>
> thank you
> Joe
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>