[R] Request for help with nls error

[Trena Phipps]

Hi.  I am attempting to run my code to get estimates for a nonlinear 
model.  Unfortunately, when I run the code, I keep getting the following 
errors (they switch back and forth depending on when I run it):

Error in nlsModel(formula, mf, start, wts) :
        singular gradient matrix at initial parameter estimates

Error in nls(ywt_vib ~ weightfunc_vib(time_vib, Wmax_star, tau_star, 
gamma_star,  :
        step factor 0.000488281 reduced below 'minFactor' of 0.000976563

I'll put my code at the bottom of this email.  Any help you can offer 
regarding this error would be greatly appreciated.  Thank you so much!!

Trena Phipps

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

vib <- matrix(scan(),ncol=3,byrow=TRUE)
1  30.0  41.1
1  45.0  51.5
1  60.0  54.0
1  75.0  56.2
1  90.0  59.1
1 120.0  59.3

time_vib <- vib[,2]
diss_vib <- vib[,3]

outfile <- "work"
cat("WORK","\n",file=outfile,"\n",append=TRUE)

n_vib <- length(time_vib)
p_vib <- 3

cmax <- 30
WLSiter <- 500
beta0_vib <- list(Wmax_star= 30, tau_star=15, gamma_star=.8)
t0 <- 20
tol <- 10^-8

meanfunc_vib <- function(time_vib,Wmax_star,tau_star,gamma_star){
   diff<-time_vib-t0
   f_vib<-Wmax_star*(1-exp(-log(2)*((diff/tau_star)^gamma_star)))
   f_vib

# meanfunc_vib_deriv <- function(time_vib,Wmax_star,tau_star,gamma_star)
meangrad_vib <- array(0,c(n_vib,p_vib))
   meangrad_vib[,1] <- 1-exp(-log(2)*(diff/tau_star)^gamma_star)
   meangrad_vib[,2] <- 
-(Wmax_star*log(2)*((diff/tau_star)^gamma_star)*exp(-log(2)*((diff/tau_star)^gamma_star)))/tau_star
   meangrad_vib[,3] <- 
Wmax_star*log(2)*((diff/tau_star)^gamma_star)*log(diff/tau_star)*exp(-log(2)*((diff/tau_star)^gamma_star))
   attr(f_vib,"gradient") <- meangrad_vib
   f_vib
}

unweightfunc_vib <- function(time_vib,Wmax_star,tau_star,gamma_star){
   diff<-time_vib-t0
   f_vib<-Wmax_star*(1-exp(-log(2)*(diff/tau_star)^gamma_star))
   f_vib
}

weightfunc_vib <- function(time_vib,Wmax_star,tau_star,gamma_star,wt_vib){
   diff<-time_vib-t0
   f_vib<-Wmax_star*(1-exp(-log(2)*(diff/tau_star)^gamma_star))
   w12_vib <- sqrt(wt_vib)
   weightf_vib <- f_vib*w12_vib
   weightgrad_vib <- 
array(0,c(n_vib,p_vib),list(NULL,c("Wmax_star","tau_star", "gamma_star")))
   weightgrad_vib[,"Wmax_star"] <- 
w12_vib*(1-exp(-log(2)*(diff/tau_star)^gamma_star))
   weightgrad_vib[,"tau_star"] <- 
-w12_vib*((Wmax_star*log(2)*((diff/tau_star)^gamma_star)*exp(-log(2)*(diff/tau_star)^gamma_star))/tau_star)
   weightgrad_vib[,"gamma_star"] <- 
w12_vib*(Wmax_star*log(2)*((diff/tau_star)^gamma_star)*log(diff/tau_star)*exp(-log(2)*(diff/tau_star)^gamma_star))
   attr(weightf_vib,"gradient") <- weightgrad_vib
   weightf_vib
}

trkfunc_vib <- function(resid_vib,mudot_vib,mu_vib,theta_vib){
   trk_vib <- resid_vib*((mudot_vib/mu_vib)**theta_vib)
   trkgrad_vib <- array(0,c(length(mu_vib),1),list(,NULL,c("theta_vib")))
   trkgrad_vib[,"theta_vib"] <- trk_vib*log(mudot_vib/mu_vib)
   attr(trk_vib,"gradient") <- trkgrad_vib   #  analytic derivative
   trk_vib
}
diss_vib_update<-diss_vib-30.50973
vib_dat<-data.frame(time_vib,diss_vib_update)

ols.fit_vib <- nls(diss_vib_update ~ 
meanfunc_vib(time_vib,Wmax_star,tau_star,gamma_star),vib_dat,beta0_vib)
bols_vib <- coef(ols.fit_vib)
mu_vib <- unweightfunc_vib(time_vib,bols_vib[1],bols_vib[2],bols_vib[3])
resid_vib <- diss_vib_update-mu_vib
sigols_vib <- sqrt(sum(resid_vib*resid_vib)/(n_vib-p_vib))

cat("FIT BY GLS",file=outfile,"\n",
     append=F)
   
cat("OLS estimate - VIB ",round(bols_vib,6),file=outfile,"\n","\n",append=T)
cat("Estimate of assumed constant SD of y 
",round(sigols_vib,6),file=outfile,
  "\n","\n","\n",append=T)

bgls_vib <- bols_vib


diss_vib_update2<-diss_vib-30.57897
vib_dat<-data.frame(time_vib,diss_vib_update2)

for (k in 1:cmax){
mu_vib <- unweightfunc_vib(time_vib,bgls_vib[1],bgls_vib[2],bgls_vib[3])
mudot_vib <- prod(mu_vib)**(1/n_vib)
mudot_vib <-rep(mudot_vib,n_vib)
resid_vib <- diss_vib_update2-mu_vib
dummy_vib <- rep(0,length(time_vib))
pldat_vib <- data.frame(resid_vib,mu_vib,mudot_vib)

theta_pl_est_vib <- nls(dummy_vib ~ 
trkfunc_vib(resid_vib,mudot_vib,mu_vib,theta_vib),pldat_vib,list(theta_vib=.8),nls.control(maxiter=300))
theta_vib <- coef(theta_pl_est_vib)

mu_vib <- unweightfunc_vib(time_vib,bgls_vib[1],bgls_vib[2],bgls_vib[3])
wt_vib <- 1/(mu_vib^(2*theta_vib))
ywt_vib <- diss_vib_update2*sqrt(wt_vib)

thedat2_vib <- data.frame(time_vib,ywt_vib,wt_vib)
       
gls_est_vib <- nls(ywt_vib ~ 
weightfunc_vib(time_vib,Wmax_star,tau_star,gamma_star,wt_vib),thedat2_vib,beta0_vib,nls.control(maxiter=200))

bgls_vib <- coef(gls_est_vib)

cat("VIB - Iteration ",k,"\n",file=outfile,append=TRUE)
cat("Estimate of theta - VIB 
",round(theta_vib,8),"\n",file=outfile,append=TRUE)
cat("GLS estimate of beta - VIB 
",round(bgls_vib,8),"\n","\n",file=outfile,append=TRUE)
}

mu_vib <- unweightfunc_vib(time_vib,bgls_vib[1],bgls_vib[2],bgls_vib[3])
resid_vib <- diss_vib_update2-mu_vib
g_vib <- mu_vib^theta_vib

sigma2_vib <- sum((resid_vib/g_vib)**2)/(n_vib-p_vib)
sigma_vib <- sqrt(sigma2_vib)

cat("Final estimate of sigma^2 - 
VIB",round(sigma2_vib,10),"\n","\n",file=outfile,append=TRUE)

sink(outfile,append=TRUE)
print(summary(gls_est_vib))
sink()