search for: varpower

...discovering that the variation in my random effect is very small. As a result, I would like to replace it as a fixed effect (i.e. essentially fit the same model but with no random effect). As I understand it I could do this using nls(), but I'm using a number of options such as weights = varPower() which I am at a loss on how to implement in that framework. Is there a way to use nlme but with out a random effect? (a bit absurd, I realize, but I have the syntax working...) Alternatively, is there a way to use "weights = varPower()" with nls? Any help would be appreciated. Sin...

Hi, I would like to use a weighted lm model to reduce heteroscendasticity. I am wondering if the only way to generate the weights in R is through the laborious process of trial and error by hand. Does anyone know if R has a function that would automatically generate the weights need for lm? Thanks, -- Tom [[alternative HTML version deleted]]

Back on January 16, a message on R-help from Ravi Varadhan described a problem with gnls using weights=varPower(). The problem was that the fit failed with error Error in eval(expr, envir, enclos) : Object "." not found I can reliably get this error in version 2.0.1-patched 2004-12-09 on Windows XP and 2.0.1-Patched 2005-01-26 on Linux. The key feature of that example is that the data are being...

I want to compare the following two models in AIC (Treat, Spotter are categorial, p is pressure, Pain is continuous) PainW.lme<-lme(Pain~p+Treat*Spotter,data=saw,random=~p|Pat, weights=varPower(form=~Pain)) # AIC= -448 Pain.lme<-lme(Pain~p+Treat*Spotter,data=saw,random=~p|Pat) #AIC = -19.7 Note the huge differences in AIC, and the estimated power of 6. A plot of the residual does not show an unusual patterns for both models. I do not trust the varPower result, but don't have any...

Hello, Excuse me for posting two questions in one day, but I figured it would be better to ask my questions in separate emails. I will again give the caveat that I'm not a statistician by training, but have a fairly decent understanding of probability and likelihood. As before, I'm trying to fit a nonlinear model to a dataset which has two main factors using nlme. Within the dataset

...ollowing structure: model<-lme(response~Covariate+TreatmentA+TreatmentB,random=~1|Block/Plot,method="ML") When I plot the residuals against the fitted values, I see a clear positive trend (meaning that the variance increases with the mean). I tried to solve this issue using weights=varPower(), but it doesn?t change the residual plot at all. How would you implement such a positive trend in the variance? I?ve tried glmmPQL (which works great with poisson errors), but using glmmPQL I can?t do model simplification. Many thanks for your help! Regards Chris.

...nd maybe this should go to R-help but here goes. I am writing a set of functions for calibration and prediction, and to calculate standard errors and intervals I need the variance function to be evaluated at new prediction points. So for instance fit<-gnls(Y~SSlogis(foo,Asym,xmid,scal),weights=varPower()) fit2<-gnls(Y~SSlogis(foo,Asym,xmid,scal),weights=varPower(form=~foo)) Now using fit or fit2 I would like to get the variance function evaluated at new points. I have played with getCovariateFormula, and looked at Initialize.gnls, summary etc. but it is not clear to me how to evaluate the fo...

...s(For3k, k=7) + pressure+ humidity, family=negbin(c(1,10)), data=efuscus) My gam.check gave me the attached result. In order to deal with my heterogeneity, I need to switch over to a gamm structure and use at least one, but possibly multiple, variance structures, and I am starting by applying varPower to my temperature covariate. (Efuscus is my square root transformed response variable). Here is the code I have for the gamm: K1 <-(efuscus~s(mic, k=7) +temp +s(date)+s(For3k, k=7) + pressure+ humidity+ windspeed + year) M17.4A <-gamm(K1, method="REML", family=negbin(c(1,10),...

...glmmPQL the Variance Structure is given as " fixed weights, Formula: ~invwt". The script shows that the function varFixed() is used, though the place where 'invwt' is defined remains unclear to me. I wonder if there is an easy way to specify another variance structure (eg varPower, etc..), preferably using an lme object of the varFunc classes ? Some trials show that the 'weights' argument of glmmPQL is just the same as in glm (which is clearly stated in the help) and I wonder actually, if not a nonsense, how to pass eg a 'weights' arguments as used in lme...

Hi, I need to correct for ar(1) behavior of my residuals of my model. I noticed that there are multiple gls models in R. I am wondering if anyone has experience in choosing between gls models. For example, how should one decide whether to use lm.gls in MASS, or gls in nlme for correcting ar(1)? Does anyone have a preference? Any advice is appreciated! Thanks, -- Tom [[alternative HTML

...mpt to use the "gnls" function and try to estimate the variance function, as a power function, I get the following error message: > ans51g <- gnls(log(b51) ~ p0 + p1/(1 + exp(-(log(dose)-p2)/p3))^p4, start=list(p0=3,p1=1,p2=4,p3=2,p4=1.5),control=gnlsControl(tol=1.e- 07),weights=varPower()) Error in eval(expr, envir, enclos) : Object "." not found > What am I doing wrong here and how can I do a GLS analysis with a variance function that is estimated from the data? Here is my data: > b51 <- c(17447.60674, 7060.37234, 2872.53012, 796.40426, 454.47222, 260....

Hello, I am trying to design new variance structures by using fixed effects variables in combination with the VarPower function. That is, I would like to create and evaluate my own variance function in the data frame and then incorporate it into the model using varPower, with value=.5. As a start, I am trying to recreate the function of VarConstPower by introducing two new variables in the data frame, d1 and d2...

Dear R members, I tried to fit an lme model and to use the glht function of multcomp. However, the glht function gives me some errors when using weights=varPower(). The glht error makes sense as glht needs factor levels and the model works fine without weights=. Does anyone know a solution so I do not have to change the lme model? Thanks Sibylle --> works fine ME$Diversity=factor(ME$Diversity) H08_lme<-lme(log(Height2005_mean)~Diversity, data=ME...

...fm1) Denom. DF: 305 numDF F-value p-value (Intercept) 1 354.7375 <.0001 sin(2 * pi * Time) 1 18.5035 <.0001 cos(2 * pi * Time) 1 1.6633 0.1981 > # variance changes with a power of the absolute fitted values? > fm2 <- update(fm1, weights = varPower()) Error in update.gls(fm1, weights = varPower()) : subscript out of bounds Execution halted (R CMD check MASS is also failing on an example that uses update.gls.) On my machine at home I run into troubles even faster when running `R CMD check nlme'. One of the early examples chokes on th...

  Hola,   Por favor, ¿podríais indicarme qué recursos (librerías o ideas) pueden resultar de utilidad para crear un gráfico del estilo del de la figura 3.8 del siguiente link?   http://www.tsc.uvigo.es/BIO/Bioing/ChrLDoc3.html#3.5   Actualmente estoy utilizando funciones muy básicas y la verdad es que no me encuentro muy satisfecha con el resultado.   Muchas gracias.   Eva [[alternative HTML

...including square variables in lme function. I've tried to work on Dialyzer data of Pinheiro and Bates'book. We fit the heteroscedastic model with: > data(Dialyzer) > fm2Dial.lme<-lme(rate~(pressure+pressure^2+pressure^3+pressure^4)*QB, + Dialyzer,~pressure+pressure^2,weights=varPower(form=~pressure)) We Obtain > fm2Dial.lme Linear mixed-effects model fit by REML Data: Dialyzer Log-restricted-likelihood: -488.4535 Fixed: rate ~ (pressure + pressure^2 + pressure^3 + pressure^4) * QB (Intercept) pressure QB300 pressure:QB300 39.362769 1....

Hi I use gnls to fit non linear models of the form y = alpha * x**beta (alpha and beta being linear functions of a 2nd regressor z i.e. alpha=a1+a2*z and beta=b1+b2*z) with variance function varPower(fitted(.)) which sounds correct for the data set I use. My purpose is to use the fitted models for predictions with other sets of regressors x, z than those used in fitting. I therefore need to estimate y with (95%) confidence intervals. Does any body knows how to do this with R ? Thanks

...ysPT,b$fit+b$se.fit*1.96,lty=2,lwd=1.5) lines(p.d$DaysPT,b$fit-b$se.fit*1.96,lty=2,lwd=1.5) points(DaysPT,Diff) However, when I add a correlation structure and/or a variance structure so that the model may look like:- gamm(Diff~s(DaysPT3)+AirToC,correlation=corCAR1(form=~DaysPT|Animal),weights=varPower(form=~DaysPT),method="REML") I get this message at the point of inputting the line "b<-predict.gam(g.m$gam,p.d,se=TRUE)" Error in model.frame(formula, rownames, variables, varnames, extras, extranames, : variable lengths differ (found for 'DaysPT') In ad...

...lain the meaning of the warning, Singular precision matrix in level -1, block 1 ? Or how to track down where it comes from? More precisely, using the nlme package, I'm issued with the warning itt2 <- lme(lrna~rx.nrti+lbrna, random=~1|patid, cor=corExp(form=~days|patid,nugget=T), weights=varPower( form=~lbrna),data=rna3) Warning messages: 1: Singular precision matrix in level -1, block 1 2: Singular precision matrix in level -1, block 1 the output is: Linear mixed-effects model fit by REML Data: rna3 Log-restricted-likelihood: -4990.142 Fixed: lrna ~ rx.nrti + lbrna (Intercept) r...

...erosced.) set.seed(1001) n=1000 x = sort(runif(n)) grow_det = exp(-2*x) grow_var = 0.1*grow_det^2 y = rnorm(n,mean=grow_det,sd=sqrt(grow_var)) dat = data.frame(x=x,y=y) ## nlme likes to have data= specified ## fit true model library(nlme) g1 = gnls(y~a*exp(-b*x), start=list(a=1,b=2), weights=varPower(form=~fitted(.)), data=dat) expdev = -2*logLik(g1) ## Fitting the true model recovers ## the true parameters nicely: coef(g1) ## Fit a series of splines: sfit <- function(d) { form <- bquote(y~ns(x,df=.(d))) gls(eval(form), weights=varPower(form=~fitted(.)),...