search for: fm1

Hi all, To follow up on an older thread, it was suggested that the following would produce confidence intervals for the estimated BLUPs from a linear mixed effect model: OrthoFem<-Orthodont[Orthodont$Sex=="Female",] fm1OrthF. <- lmer(distance~age+(age|Subject), data=OrthoFem) fm1.s <- coef(fm1OrthF.)$Subject fm1.s.var <- fm1OrthF. at bVar$Subject fm1.s0.s <- sqrt(fm1.s.var[1,1,]) fm1.s0.a <- sqrt(fm1.s.var[2,2,]) fm1.s[,1]+outer(fm1.s0.s, c(-2,0,2)) fm1.s[,2]+outer(fm1.s0.a, c(-2,0,2)) [,1...

...but I'm not sure I understood. Is the approach shown below correct? Rick B. # Orthodont is from nlme (can't have both nlme and lme4 loaded at same time!) # OrthoFem<-Orthodont[Orthodont$Sex=="Female",] # http://tolstoy.newcastle.edu.au/R/help/06/03/23758.html library(lme4) fm1OrthF. <- lmer(distance~age+(age|Subject), data=OrthoFem) lmer(distance~age+(age|Subject), data=OrthoFem)@bVar$Subject[2,2,]* (attr(VarCorr(lmer(distance~age+(age| Subject),data=OrthoFem)),"sc")^2)[1] (attr(VarCorr(fm1OrthF.),"sc")^2)[1] fm1.s <- coef(fm1OrthF.)$...

...specific question in r-help etiquette? # example set.seed(123) x <- arima.sim(list(order = c(1,0,0), ar = 0.7), n = 100) y <-x + arima.sim(list(order = c(1,0,0), ar = 0.7), n = 100) x <- c(x) y <- c(y) lm1 <- lm(y~x) ar(residuals(lm1)) # indicates an ar1 model cs1 <- corARMA(p=1) fm1 <- gls(y~x,corr=cs1) summary(fm1) # get fits fits <- predict(fm1) # use coef to get fits fits2 <- coef(fm1)[1] + (coef(fm1)[2] * x) plot(fits,fits2)

Hi all, I am having issues comparing models with lmer. As an example, when I run the code below the model summaries (AIC, BIC, loglik) differ between the summary() and anova() commands. Can anyone clear up what's wrong? Thank you! Darren Ward library(lme4) data(sleepstudy) fm1<-lmer(Reaction ~ Days + (1|Subject), sleepstudy) summary(fm1) fm2<-lmer(Reaction ~ Days + (Days|Subject), sleepstudy) summary(fm2) anova(fm1, fm2) Sample output: > summary(fm1) Linear mixed-effects model fit by REML Formula: Reaction ~ Days + (1 | Subject) Data: sleepstudy AIC B...

Hi List, I have to compare a relationship between y and x for two locations. I found logistic regression fits both datasets well, but I am not sure how to test if relationships for both sites are significantly different. I searched the r site, however no answers exactly match the question. I used Tukey's HSD to compare two means, but the relationship in my study was not simply linear. So I

Hello, I am looking for a way to obtain standard errors for emprirical Bayes estimates of a model fitted with lmer (like the ones plotted on page 14 of the document available at http://www.eric.ed.gov/ERICDocs/data/ericdocs2/content_storage_01/0000000b/80/2b/b3/94.pdf). Harold Doran mentioned (http://tolstoy.newcastle.edu.au/~rking/R/help/05/08/10638.html) that the posterior modes' variances

...in Australia) is as follows: CONDT <- c(0, 1, 1, 0, 2, 3, 0, 1, 1, 1, 1, 2, 0, 1, 3, 0, 1, 2, 1, 3, 3, 4, 1, 3, 2, 0, 2, 0, 3, 0, 0, 1, 1, 1, 1, 0, 0, 2, 2, 0, 1, 2, 0, 0, 1, 1, 1, 0, 2) The Poisson model estimating the mean: summary(fm1 <- glm(CONDT ~ 1, family = poisson(link = 'identity'))) The estimated coefficient obtained matches the results in Dobson exactly. But I am unable to replicate the value for the log-likelihood using the results on deviance or null deviance (the same here) from R. If you plug into the log...

...4374072/test_text.txt test_text.txt (see upload) CODE: ( which occurs after TEST <- read.table("test_text.txt") ) i <- 2; sink (file = "test_output.txt", append = FALSE) mydf <- data.frame(TEST) for (j in 1:ncol(mydf)-1) { var1 <- mydf[,1] var2 <- mydf[,i] fm1 <- aov(var1 ~ var2) tky <- TukeyHSD(fm1) otpt <- capture.output(summary(fm1)) i <- i+1; lines <- as.vector(unlist(strsplit(otpt[2]," ")),mode="list") # gets the p-value if (grepl("[1234567890]",lines[14],perl = TRUE)) { #make sure that slot for p-va...

...ummary(lm) help. There seems to be a discrepancy between the what R produced and the manual calculation. Does anyone know why this is so? What does the multiple r squared reported in summary(lm) represent? # The test case: x <- c(1,2,3,4) y <- c(1.6,4.4,5.5,8.3) dummy <- data.frame(x, y) fm1 <- lm(y ~ x-1, data = dummy) summary(fm1) betax <- fm1$coeff[x] * sd(x) / sd(y) # cd is coefficient of determination cd <- betax * cor(y, x) Thanks.

...em. If you have this problem, you can either exit Rcmdr before calling TukeyHSD or you can set na.action to na.omit. The code below demonstrates the situation. Cheers, Bob data(warpbreaks) head(warpbreaks) # Introduce a missing value: warpbreaks$breaks[1] <- NA head(warpbreaks) # Do a model: fm1 <- aov(breaks ~ tension, data = warpbreaks) TukeyHSD(fm1, "tension", ordered = TRUE) # Setting na.exclude or starting Rcmdr will kill the confidence intervals: options(na.action = na.exclude) fm1 <- aov(breaks ~ tension, data = warpbreaks) TukeyHSD(fm1, "tension", ordered...

...e levels of Subj indicate a grouping structure (k subjects) and Trt is a two-level factor (two treatments) for which there are several (n) responses y from each treatment and subject combination. If one suspects a subject by treatment interaction, either of the following models seem natural > fm1 <- lme(y ~ Trt, random = list(Subj = pdDiag(~ Trt))) > fm2 <- lme(y ~ trt, random = ~ 1 | Subj/Trt) These models seem to correspond to the same situation. Both have two variance components (subject and treatment within subject). However, they result in different denominator degrees of fre...

Dear R-Users, Is it possible to get the covariance matrix from an lme model that did not converge ? I am doing a simulation which entails fitting linear mixed models, using a "for loop". Within each loop, i generate a new data set and analyze it using a mixed model. The loop stops When the "lme function" does not converge for a simulated dataset. I want to

...ned using nlme function gls. The anova procedure fails to locate the second of the objects. The following code, borrowed from the help page of anova.gls, exemplifies: --------------- start example code --------------- library(nlme) ## stolen from example(anova.gls) # AR(1) errors within each Mare fm1 <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), Ovary, correlation = corAR1(form = ~ 1 | Mare)) anova(fm1) # variance changes with a power of the absolute fitted values? fm2 <- update(fm1, weights = varPower()) anova(fm1, fm2) ## now define a little function dummy <- functi...

...n see my attempts together with results. I hope somebody can point me to right direction. I am probably somewhere close but I have no clue, which parameter I shall modify to get measured points, fitted lines and vertical lines in panels together. Please help Thank you Best regards. Petr Pikal fm1 <- lme(Orthodont) # standard plot plot(augPred(fm1, level = 0:1, length.out = 2)) #plot with vertical but without points and fitted lines plot(augPred(fm1, level = 0:1, length.out = 2), panel=function(v,...) { panel.abline(v=10)} ) # plot with vertical but without fitted lines plot(augPred(fm...

...ing to understand the model specification formalities for 'lme', and the documentation is leaving me a bit confused. Specifically, using the example dataset 'Orthodont' in the 'nlme' package, first I use the invocation given in the example shown by "?lme": > fm1 <- lme(distance ~ age, data = Orthodont) # random is ~ age Despite the Comment ("# random is ~ age"), > summary(fm1) says that [...] Random effects: Formula: ~age | Subject Structure: General positive-definite [...] Fixed effects: distance ~ age In view of the...

When I try to run the example from Variogram with an lme object, I get an error (although summary works): R : Copyright 2005, The R Foundation for Statistical Computing Version 2.2.1 (2005-12-20 r36812) ISBN 3-900051-07-0 ... > fm1 <- lme(weight ~ Time * Diet, BodyWeight, ~ Time | Rat) Error: couldn't find function "lme" > Variogram(fm1, form = ~ Time | Rat, nint = 10, robust = TRUE) Error: couldn't find function "Variogram" > library(nlme) > fm1 <- lme(weight ~ Time * Diet, Bo...

Dear all, I have three concerns: 1) I am running models with the lme4 package. I cannot find a way to pull out a vector of the fitted values and the residuals. Does anybody know how to do it? 2) How can I nest a random effect variable into a "two-level" fixed effect variable? 3) Suppose I have the following model: y = a + b|c + d + error, where 'a' is a fixed effect, 'c'

...ail Rail travel 1 1 55 2 1 53 3 1 54 4 2 26 5 2 37 6 2 32 7 3 78 8 3 91 9 3 85 10 4 92 11 4 100 12 4 96 13 5 49 14 5 51 15 5 50 16 6 80 17 6 85 18 6 83 > > fm1 <- lme(travel~1, data=Rail, random=~1|Rail, method="ML") > summary(fm1) Linear mixed-effects model fit by maximum likelihood Data: Rail AIC BIC logLik 134.5600 137.2312 -64.28002 Random effects: Formula: ~1 | Rail (Intercept) Residual StdDev: 22.62435...

I am having a problem extracting from "mer" objects.    I have constructed my problem using existing datasets.   Using the following commands:   require(lme4) fm1 <- lmer(Yield ~ 1 + (1 | Batch), Dyestuff) fixef(fm1) I get the following error message: "Error in UseMethod("fixef") : no applicable method for "fixef""   I know that "fixef" is in lme4 and UseMethod is in "base", so I am stuck trying to sort o...

dear all, I would like to get the lme call without fitting the relevant model. library(nlme) data(Orthodont) fm1 <- lme(distance ~ age, random=list(Subject=~age),data = Orthodont) To get fm1$call without fitting the model I use call(): my.cc<-call("lme.formula", fixed= distance ~ age, random = list(Subject = ~age)) However the two calls are not the same (apart from the data argument I am...