search for: modelstruct

...,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 ap...

The print method for lme *prints out* the fitted correlation matrix for the random effects. Is there any way to get these values as an object in R? I have examined the components of the lme object (called "junk" in the example below) and the components of summary(junk) without finding these numbers. (How I did this: I dumped the entire lme object to a text file and then used egrep to

...nction, but sometimes when the estimated Hessian matrix is not positive definite or something like that (i am not quite sure), the intervals function delivers an error message. Thus, i will like to ask if there is another way to access these values. I tried using the following code: hh$modelStruct$corStruct[1] hh$modelStruct$varStruct[1] Rather the output was: > hh$modelStruct$corStruct[1] [1] -1.308580 > hh$modelStruct$varStruct[1] [1] -0.03124255 I presume there is a way to calculate the correlation and variance function coefficients using these values. Cou...

...nction, but sometimes when the estimated Hessian matrix is not positive definite or something like that (i am not quite sure), the intervals function delivers an error message. Thus, i will like to ask if there is another way to access these values. I tried using the following code: hh$modelStruct$corStruct[1] hh$modelStruct$varStruct[1] Rather the output was: > hh$modelStruct$corStruct[1] [1] -1.308580 > hh$modelStruct$varStruct[1] [1] -0.03124255 I presume there is a way to calculate the correlation and variance function coefficients using these values. Cou...

....graves at pdf.com] Sent: Sunday, January 15, 2006 6:41 PM To: Rand, Hugh Cc: 'r-help at lists.R-project.org' Subject: Re: [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...

Dear all, Here is the question: For example, using the "petrol" data offered with R. pet3.lme<-lme(Y~SG+VP+V10+EP,random=~1|No,data=petrol) pet3.lme$sigma gives the residual StdDev. But I can't figure out how to extract the "(intercept) StdDev", although it is in the print out if I do "summary(pet3.lme)". In

...tiple imputation for missing values of the response variable in the model. I am trying to extract the within group correlation matrix or covariance matrix. here is my code: f = lme(y ~x0+x1+trt+tim+x1:tim +tim:trt,random=~-1|subj, data=dat,corr =corAR1()) > f$sigma [1] 2.330854 b=summary(f)$modelStruct$corStruct > b Correlation structure of class corAR1 representing Phi 0.8518711 I think 0.8518711 and f$sigma is what I need to reconstruct the variance-cavariance matrix. So, How can I extract 0.8518711 so that I can assign it to a variable? Also, I don't understand what the follow...

...fit my simulated data, we use varConstPower ( constant plus power variance function). fm<-lme(UPDRS~time,data=data.simula,random=~time,method="ML",weights=varConstPower(fixed=list(power=1))) I extract the results of this function by using the following codes: y<-summary(fm) x<-y$modelStruct$varStruct > x Variance function structure of class varConstPower representing const power 1.184535e-15 1.000000e+00 These are the constants and power that apppear in the summary(fm), and that I want to extract. Hower I got different value of const when extracting from x: >...

Hi, lme() users, Can some one tell me how to do this. I model Orthodont with the same G for random variables, but different R{i}'s for boys and girls, so that I can get sigma1_square_hat for boys and sigma2_square_hat for girls. The model is Y{i}=X{i}beta + Z{i}b + e{i} b ~ iid N(0,G) and e{i} ~ iid N(0,R{i}) i=1,2 orth.lme <- lme(distance ~ Sex * age, data=Orthodont, random=~age|Subject,

...servations: 291 Number of Groups: 119 > When I try to get the contrast I get one of two errors each time. Trying for a paired contrast... tst <- contrast(lme.fnl.REML, a=list(Yr=levels(Yr), Tr="R"), b=list (Yr=levels(Yr, Tr="T97"))) Error in gendata.default(fit = list(modelStruct = list(reStruct = list ( : not enough factors Trying to include the other factor to make R happy.... > tst <- contrast(lme.fnl.REML, a=list(Yr=levels(Yr), Tr="R"), b=list (Yr=levels(Yr, Tr="T97")), c=list(Yr=levels(Yr, Tr="T98")) + ) Error in contrastCalc(f...

...=sqrt(0.5) alpha0=0.1 alpha1=3 m=200 n=200 c0=-2.413; nu0=-0.2; nu1=0.3 simu.dat=simulate.pilot(m, mn, sigma, alpha0, alpha1, c0, nu0, nu1) fit1=try(gls(Y~W, data=simu.dat, weights=varComb(varExp(form=~W), varExp(form=~I(W^2))))) c0.hat= log(fit1$sigma^2) nu0.hat=2*fit1$modelStruct$varStruct$A[1] nu1.hat=2*fit1$modelStruct$varStruct$B[1] > c0.hat [1] -1.570104 > nu0.hat [1] -0.787264 > nu1.hat [1] 0.4057129 > Thanks Xuesong CONFIDENTIALITY NOTICE: This e-mail message, including a...{{dropped:12}}

...} mn=3.3 sigma=sqrt(0.5) alpha0=0.1 alpha1=3 m=200 n=200 c0=-2.413; nu0=-0.2; nu1=0.3 simu.dat=simulate.pilot(m, mn, sigma, alpha0, alpha1, c0, nu0, nu1) fit1=try(gls(Y~W, data=simu.dat, weights=varComb(varExp(form=~W), varExp(form=~I(W^2))))) c0.hat= log(fit1$sigma^2) nu0.hat=2*fit1$modelStruct$varStruct$A[1] nu1.hat=2*fit1$modelStruct$varStruct$B[1] > c0.hat [1] -1.570104 > nu0.hat [1] -0.787264 > nu1.hat [1] 0.4057129 > Thanks Xuesong

...res, but I thought that maybe optimum weights from gls might be used as weights in lm (as shown below), but apparently this is not the case. See: library(nlme) f1 <- gls(Petal.Width ~ Species / Petal.Length, data = iris, weights = varIdent(form = ~ 1 | Species)) aa <- attributes(summary(f1)$modelStruct$varStruct)$weights f2 <- lm(Petal.Width ~ Species / Petal.Length, data = iris, weights = aa) summary(f1)$tTable; summary(f2) So, the two models with the very same weights do differ (in terms of standard errors). Could you please explain why? Are these different types of weights? Many thanks i...

...e are presented by print.lme as: ... (Intercept) Residual StdDev: 1.410635 0.7800512 ... (data taken from ?lme's examples section) I can get the residuals with x$sigma where x is the fitted lme object. I can't find the intercept, though. The closest I've got is: > tmp$modelStruct$reStruct Random effects: Formula: ~1 | ID (Intercept) Residual StdDev: 1.808388 1 which still doesn't give me just the number. I've been looking in the summary() and print() methods to the structures, but haven't had any luck in pinpointing where this number is...

...ach level's equation, and residual is the residual variance. The problem is getting access to var1..varn (or, as R generally reports it, sd1..sdn) programmatically. I can get the residual standard deviation with: print(model$sigma) and I can view the standard deviations with: summary(model$modelStruct$lmeStruct) But I can't figure out how to get at the standard deviations without the human-readable print format. Thanks- Andy Perrin ---------------------------------------------------------------------- Andrew J Perrin - http://www.unc.edu/~aperrin Assistant Professor of Sociology, U of Nor...

...ormula, and looked at Initialize.gnls, summary etc. but it is not clear to me how to evaluate the form component, especially in the case of fit above where form=~fitted(.), in any safe way. I can grep for "fitted" in the formula eg. grep("fitted",deparse(getCovariateFormula(fit$modelStruct$varFunc))) and try to calculate predicted values from the model for the new points but how to substitute back in the new terms? I don't need this problem solved on a platter, I just need to unedrstand an approach, because my stabs are failing. Thanks Nicholas

Greetings, 1) Is there a nice way of extracting the variance estimates from an lme fit? They don't seem to be part of the lme object. 2) In a series of simulations, I am finding that with ML fitting one of my random effect variances is sometimes being estimated as essentially zero with massive CI instead of the finite value it should have, whilst using REML I get the expected value. I guess

...and reasonable to try and get the: - null deviance - the total deviance due to the the random effect(s) - the residual deviance? With the idea that Null deviance = Fixed effects + Random Effects + Residuals If yes how to do it ? A lme object provides the following: > names(glm6) [1] "modelStruct" "dims" "contrasts" "coefficients" [5] "varFix" "sigma" "apVar" "logLik" [9] "numIter" "groups" "call" "method" [13] &quo...

...do? The simplest >thing, if the method is giving incorrect results, is to eliminate the >method entirely. Any objections to my doing that? > >The method is currently defined as > >getVarCov.gls <- > function(obj, individual = 1, ...) >{ > S <- corMatrix(obj$modelStruct$corStruct)[[individual]] > if (!is.null( obj$modelStruct$varStruct)) > { > ind <- obj$groups==individual > vw <- 1/varWeights(obj$modelStruct$varStruct)[ind] > } > else vw <- rep(1,nrow(S)) > vars <- (obj$sigma^2)*vw > resu...

...eatment among forests levels) I used function contrast (package contrast) whit following script to probe the hypotesis 1.: con<-contrast(FINAL, list(Lugar= 'Xmatkuil', Tra=1), list(Lugar='Xmatkuil', Tra = 0)) but i found this error message: Error en gendata.default(fit = list(modelStruct = list(varStruct = c(-0.253689933940456, : not enough factors I would be grateful if somebody tell me I'm doing wrong with my contrast function script. Thanks in advance, Marylin Bejarano PHd candidate in Ecology Institute of Mexico's National Autonomous University [[alternative HT...