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()