[R] Problems and bugs in vgam()

[John Poulsen]

Hello R-Users,

I have recently run into several problems using vgam() in the VGAM 
package.  I am hoping someone might have some solutions...

Briefly, I have been trying to fit GAM models for zero-altered negative 
binomial models.

1. When fitting smoothed parameters (e.g. s(X, df=2)) changing the 
degrees-of-freedom has no effect on the level of smoothing (e.g. number 
of knots for the spline).  This appears to be the case even for the 
example on the vgam help page.

data(hunua)
fit2 = vgam(agaaus ~ s(altitude, df=2), binomialff, hunua)

 > coef(fit2)
         (Intercept) s(altitude, df = 2)
       -1.1661259280        0.0003932463

fit10 = vgam(agaaus ~ s(altitude, df=10), binomialff, hunua)

 > coef(fit10)
          (Intercept) s(altitude, df = 10)
        -1.1661259280         0.0003932463

2. There may be a bug when trying to run gam (mgcv) and vgam(VGAM) 
together.  When I run vgam first, I am unable to get smoothed parameters 
using gam.  This presents a problem if one would like to use the 
parameter estimates from a gam model to try to fit a vgam model.

3. For my own data, I have been having a very difficult time getting 
vgam to converge.  The help pages says that convergence may depend 
strongly on the initial values (coefstart) provided, and I have been 
modifying them but to no avail.  Is there another way or more values 
that need to be given to make vgam models fit?

I frequently run into the error below.  Any advice on how to overcome it?

Applying Greenstadt modification to 357 matrices
Error in ans[, index] <- tmp777 :
   number of items to replace is not a multiple of replacement length
In addition: Warning message:
In checkwz(wz, M = M, trace = trace, wzeps = control$wzepsilon) :
   357 elements replaced by 1.819e-12

In case someone has the interest or time in helping me get models to 
converge, I provide an example model and data.  (Data are attached and 
the model code is below.)  Thanks for your help.

library(VGAM)
Cnictdat=read.csv("Cnictdat.csv", header=T, sep=",")

coefnict=c(-6.249e+00, -8.049e+00, -1.0e+01, -1.66e+00,
            6.620960e-01, -1.046e-02, 4.586179e-01, -1.189e+00,
           -2.712e-03, -2.126218e+00, 7.242e-01, -1.149e+00,
            9.377e-01, -8.457e-03, -3.035807e-02, -1.697e-02,
            2.279305e-02, -4.7e-07, -1.797e-07, -1.183e-07, -4.3494e-05)
abunddst=c(Cnictdat$SegDistKm)*1000
esw=0.5

fit4=vgam(Cnict~factor(Hab)+s(Clr, spar=2)+s(Rd, spar=2)+s(YL, 
spar=2)+ 		      s(Pop, spar=2)+offset(log(2*esw*abunddst)),
               zanegbinomial, data=Cnictdat, coefstart=coefnict,
               control=vgam.control(se.fit=T, exp=1e-4, maxit=500,
	      trace=T))

Thanks,
John