R help - Jun 2014 - How to define a "nonlin.function" object for the spherical function to be used in the gnm model?

Hi,


I encountered with some difficulty in preparing the nonlin.function for the
gnm function. Please help me.


 I want to use the gnm function in the "gnm" package to fit the
following
spherical function


y =  a + b (1.5*x - 0.5 *(x/c))         if  x< c

y =  a + b                                            if  x >= c


I understand that the fit.variogram function in the gstat function was
designed for fitting variogram, but it seems to be questionable because of
the normality assumption on the residuals. It is interested that whether
the fitted model would be changed as another distribution was assumed.


Seems to me, the gnm function in the gnm/R package could be used for the
this. However, I failed to customize a proper "nonlin" object for the
aforementioned formula after a few hours of try.

The following are my crappy code defining the nonlin object.

-------------------------------------------------------------

myfun<-function(x) {

  list( predictors = list(nuggest = .1, psill = 1, range=5),    variables
list(substitute(x)),

    term <- function(predLabels, varLabels) {

      paste0("pred=", predLabels[1], "+", predLabels[2],
"*(1.0 - exp(-",
varLabels[1], "/", predLabels[3], ")")

#       paste0("pred[which(", varLabels[1], ">",
predLabels[3], ")]=",
predLabels[3], "+", predLabels[2])

#       paste0("pred")

    }

  )

}

class(myfun) <- "nonlin"

-------------------------------------------------------------

In addition, I have also tried the gnls function in the nlme package but
did not find any parameters in function would allow to specify a
distribution for the residuals.


I appreciate any help, suggestion, comment on my question.Thank you very
much.


Best regards,


Zhiqiu



The following code are largely irrelevant to my question, but it could be
used to generate a semivarogram table for the analysis.

#The Sph model in SP analysis

sph.fun<-function(dist, nugget, psill, range) {

  res<-nugget+psill*(1.5*dist/range-.5*(dist/range)^3)

  id.far=which(dist>=range)

  if (NROW(id.far)>0)  res[id.far]= sum(nugget, psill)

  res

}



#simulate semivariogram data

library(gstat)

dist=1:100

dat=data.frame(np=300.0, dist=(dist)+0.001,

               gamma=sph.fun(dist, nugget = 1, psill = 9, range
20)+rnorm(NROW(dist))*1e-3,

               dir.hor=0.0, dir.ver=0.0,  id=as.factor(1))

v=dat;

class(v)<-"gstatVariogram"



#By defalut, WLS is used and the weight is  w=np/dist^2

ft1=fit.variogram(v, vgm(nugget = 1, mod="Sph", psill = 1, range = 5))



#use the nls to fit the data

ft2=nls(gamma ~ sph.fun(dist, nugget, psill, range),

        data = dat,

        start = list(nugget=1,  psill= 1, range = 5),

        weight=np/dist^2,

        algorithm="port",

        lower=c(0, 0, 0),

        trace=F)

	[[alternative HTML version deleted]]