StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

POCannot update my nlme model using varPower or other varClasses
text
Body
copied!<p>I have fit a model using nlme and I don't think the residuals look good enough, so I am trying a transformation by</p> <pre><code>fit.nlme6 <- update(fit.nlme5, weights = varPower()) </code></pre> <p>I get <code>Error: Singularity in backsolve at level 0, block 1</code>. I have tried other classes that I don't think make much sense anyway and tried various <code>forms = ~</code> and <code>fixed</code>. All with the same error message, sometimes with bonus warning messages.</p> <p>Here are the residuals. I think <code>varPower</code> should work perfectly, so why doesn't it? </p> <p>residuals:</p> <p><img src="https://i.imgur.com/GQRgk.png"></p> <p>More information, <code>fit.nlme5</code> is a model fit to a beta growth function that has three parameters, <code>w.max</code> (maximum biomass), <code>t.e</code> (moment growth ends), and <code>t.m</code> (moment of maximum growth). The model looks like this</p> <pre><code>fit.nlme5 <- update(fit.nlme4, fixed = list(t.e ~ trt + ground + trt:ground, w.max + t.m ~ 1), start = c(cfsTD[1:4], rep(0,2) ,cfsTD[5:6])) </code></pre> <p>There are three treatments (trt) in two locations (ground). After the residuals get fixed, I'll run some contrasts to compare the treatments in the different locations.</p> <p>Here is the data <a href="https://dl.dropbox.com/u/21080842/cobsgddv8.txt" rel="nofollow noreferrer">https://dl.dropbox.com/u/21080842/cobsgddv8.txt</a></p> <p>And the code one would need to get to the final model:</p> <pre><code> #You'll need these functions ## Implementing the beta growth function from (Yin et al 2003) bgfInit <- function(mCall, LHS, data){ xy <- sortedXyData(mCall[["time"]], LHS, data) if(nrow(xy) < 4){ stop("Too few distinct input values to fit a bgf") } w.max <- max(xy[,"y"]) t.e <- NLSstClosestX(xy, w.max) t.m <- NLSstClosestX(xy, w.max/2) value <- c(w.max, t.e, t.m) names(value) <- mCall[c("w.max","t.e","t.m")] value } bgf <- function(time, w.max, t.e, t.m){ .expr1 <- t.e / (t.e - t.m) .expr2 <- (time/t.e)^.expr1 .expr3 <- (1 + (t.e - time)/(t.e - t.m)) .value <- w.max * .expr3 * .expr2 ## Derivative with respect to t.e .exp1 <- ((time/t.e)^(t.e/(t.e - t.m))) * ((t.e-time)/(t.e-t.m) + 1) .exp2 <- (log(time/t.e)*((1/(t.e-t.m) - (t.e/(t.e-t.m)^2) - (1/(t.e - t.m)))))*w.max .exp3 <- (time/t.e)^(t.e/(t.e-t.m)) .exp4 <- w.max * ((1/(t.e-t.m)) - ((t.e - time)/(t.e-t.m)^2)) .exp5 <- .exp1 * .exp2 + .exp3 * .exp4 ## Derivative with respect to t.m .ex1 <- t.e * (time/t.e)^((t.e/(t.e - t.m))) * log(time/t.e) * ((t.e - time)/(t.e - t.m) + 1) * w.max .ex2 <- (t.e - time) * w.max * (time/t.e)^(t.e/(t.e-t.m)) .ex3 <- (t.e - t.m)^2 .ex4 <- .ex1 / .ex3 + .ex2 / .ex3 .actualArgs <- as.list(match.call()[c("w.max", "t.e", "t.m")]) ## Gradient if (all(unlist(lapply(.actualArgs, is.name)))) { .grad <- array(0, c(length(.value), 3L), list(NULL, c("w.max", "t.e", "t.m"))) .grad[, "w.max"] <- .expr3 * .expr2 .grad[, "t.e"] <- .exp5 .grad[, "t.m"] <- .ex4 dimnames(.grad) <- list(NULL, .actualArgs) attr(.value, "gradient") <- .grad } .value } SSbgf <- selfStart(bgf, initial = bgfInit, c("w.max", "t.e", "t.m")) #Now for the data and fitting grow<-read.table("cobsgddv8.txt", header=T) library(nlme) grow10<-subset(grow, grow$year == "2010") grow10$EU<- with(grow10, factor(ground):factor(plot)) grow10G<-groupedData(mass ~ gdd | EU, data=grow10) fit.beta.10 <- nlsList(mass ~ SSbgf(gdd, w.max, t.e, t.m), data = grow10G) fit.nlme.10<-nlme(fit.beta.10, random=pdDiag(w.max ~1)) cfs <- fixef(fit.nlme.10) fit.nlme2 <- update(fit.nlme.10, fixed = list(t.e ~ trt, w.max + t.m ~ 1), start = c(cfs[2], rep(0,2), cfs[c(1,3)])) cfsT <- fixef(fit.nlme2) fit.nlme3 <- update(fit.nlme.10, fixed = list(t.e ~ ground, w.max + t.m ~ 1), start = c(cfs[2], rep(0,1), cfs[c(1,3)])) cfsG <- fixef(fit.nlme3) fit.nlme4 <- update(fit.nlme.10, fixed = list(t.e ~ trt + ground, w.max + t.m ~ 1), start = c(cfsT[1:2], cfsG[1:2], cfs[c(1,3)])) cfsTD <- fixef(fit.nlme4) fit.nlme5 <- update(fit.nlme4, fixed = list(t.e ~ trt + ground + trt:ground, w.max + t.m ~ 1), start = c(cfsTD[1:4], rep(0,2) ,cfsTD[5:6])) fit.nlme6 <- update(fit.nlme5, weights = varPower()) </code></pre>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload