Rand, Hugh
2006-Jan-09 15:37 UTC
[R] trouble with extraction/interpretation of variance structure para meters from a model built using gnls and varConstPower
I have been using gnls with the weights argument (and varConstPower) to specify a variance structure for curve fits. In attempting to extract the parameters for the variance model I am seeing results I don't understand. When I simply display the model (or use "summary" on the model), I get what seem like reasonable values for both "power" and "const". When I actually try to extract the values, I get the same number for the "power", but a different (and less sensible) value for "const". The simplest example I can come up with that shows the problem is as follows: #Set up data x = c(0,0,0,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7); y0 = c(.1,.1,.1, .5,.5,.5, 1,1,1, 2,2,2, 4,4,4, 7,7,7, 9,9,9, 10,10,10); yp = c(0,.03,.05, 0,.05,.01, 0,.07,.03, .1,0,.2, .2,.1,0, .3,0,.1, 0,.3,.4, .3,.5,0); y = y0 + 4*yp; data = data.frame(x=x,y=y); #Run model library(nlme) model3 = try(gnls(y ~ SSfpl(x,A,B,xmid,scal),data=data,weights=varConstPower(const=1,power=0,form~fitted(.)))); #Examine results model3; #const = .6551298, power .8913665 summary(model3); #const = .6551298, power .8913665 coef(model3$modelStruct$varStruct)["power"]; # power .8913665 coef(model3$modelStruct$varStruct)["const"]; #const = -0.4229219 coef.varFunc(model3$modelStruct$varStruct); #R can't seem to find this function, Splus can Any advice on what I am doing wrong would be appreciated. Hugh Rand Senior Scientist Amgen rand at amgen.com
Spencer Graves
2006-Jan-16 02:40 UTC
[R] trouble with extraction/interpretation of variance structure para meters from a model built using gnls and varConstPower
How about this: > exp(coef(model3$modelStruct$varStruct)["const"]) const 0.6551298 Does that answer the question about not understanding the connection between summary(model3) and coef(model3$modelStruct$varStruct)["const"]? Regarding the question about R not being able to find 'coef.varFunc', I tried the following: > methods("coef") <snip> coef.varFunc* [37] coef.varIdent* coef.varPower* Non-visible functions are asterisked Since "coef.varFunc" is "asterisked", I tried 'getAnywhere(coef.varFunc)'. I walked through this code line by line, until I found the following: > (val <- as.vector(model3$modelStruct$varStruct)) Variance function structure of class varConstPower representing const power 0.6551298 0.8913665 Answer the questions? spencer graves Rand, Hugh wrote:> I have been using gnls with the weights argument (and varConstPower) to > specify a variance structure for curve fits. In attempting to extract the > parameters for the variance model I am seeing results I don't understand. > When I simply display the model (or use "summary" on the model), I get what > seem like reasonable values for both "power" and "const". When I actually > try to extract the values, I get the same number for the "power", but a > different (and less sensible) value for "const". > > The simplest example I can come up with that shows the problem is as > follows: > > > #Set up data > > x = c(0,0,0,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7); > y0 = c(.1,.1,.1, .5,.5,.5, 1,1,1, 2,2,2, 4,4,4, 7,7,7, 9,9,9, 10,10,10); > yp = c(0,.03,.05, 0,.05,.01, 0,.07,.03, .1,0,.2, .2,.1,0, .3,0,.1, > 0,.3,.4, .3,.5,0); > y = y0 + 4*yp; > data = data.frame(x=x,y=y); > > #Run model > > library(nlme) > model3 = try(gnls(y ~ > SSfpl(x,A,B,xmid,scal),data=data,weights=varConstPower(const=1,power=0,form> ~fitted(.)))); > > #Examine results > > model3; #const = .6551298, power > .8913665 > summary(model3); #const = .6551298, power > .8913665 > > coef(model3$modelStruct$varStruct)["power"]; # power > .8913665 > coef(model3$modelStruct$varStruct)["const"]; #const = -0.4229219 > coef.varFunc(model3$modelStruct$varStruct); #R can't seem to find this > function, Splus can > > > Any advice on what I am doing wrong would be appreciated. > > Hugh Rand > Senior Scientist > Amgen > rand at amgen.com > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html