search for: fit0

...rray(runif(720),c(240,3),dimnames=list(NULL,c('x1','x2','y' )))) id=rep(1:120,2); datx=cbind(id,datx) #give x1 a slight relation with y (only necessary to make the random effects non-zero in this artificial example) datx$x1=(datx$y*0.1)+datx$x1 library(lme4) #fit the data fit0=lmer(y~x1+x2+(1|id), data=datx); print(summary(fit0),corr=F) fit1=lmer(y~x1+x2+(1+x1|id), data=datx); print(summary(fit1),corr=F) #compare the models anova(fit0,fit1) Now, look at the output, below. You can see that the AIC from "print(summary(fit0))" is 87.34, but the AIC for fit0 in...

...on the three models. I expect that the p.value for the comparison of model 1 and model 2 would remain the same, whether or not I add a third model to be compared with. However, the P values change, and I do not understand why. Here is an example code (following with it's input): data(barro) fit0 <- rq(y.net ~ lgdp2 + fse2 , data = barro) fit1 <- rq(y.net ~ lgdp2 + fse2 + gedy2 , data = barro) fit2 <- rq(y.net ~ lgdp2 + fse2 + gedy2 + Iy2 , data = barro) anova(fit0,fit1,fit2, R = 1000) anova(fit0,fit1, R = 1000) Output: > data(barro) > fit0 <- rq(y.net ~ lgdp2 + fse...

...l exist). Sample code: > ## Avoid printing to unwarranted accuracy > od <- options(digits = 5) > x <- 0:10 > y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8) > > ## Easy one-dimensional MLE: > nLL <- function(lambda, y) -sum(stats::dpois(y, lambda, log = TRUE)) > fit0 <- mle(nLL, start = list(lambda = 5), fixed=list(y=y), nobs = NROW(y)) > confint(fit0) Profiling... 2.5 % 97.5 % 9.6524 13.6716 > rm(y) > confint(fit0) Profiling... Error in eval(expr, p) : object 'y' not found In this sample code, I?m showing that after the removal of y...

I am trying to find methods for testing and visualizing calibration to Cox models with time-depended coefficients. I have read this nice article <http://journals.sagepub.com/doi/10.1177/0962280213497434>. In this paper, we can fit three models: fit0 <- coxph(Surv(futime, status) ~ x1 + x2 + x3, data = data0) p <- log(predict(fit0, newdata = data1, type = "expected")) lp <- predict(fit0, newdata = data1, type = "lp") logbase <- p - lp fit1 <- glm(y ~ offset(p), family = poisson, data = data1) fit2 <- glm(y...

Hi , I want to apply the following code fo my data with 400 predictors. I was wondering if there ia an alternative way instead of typing 400 predictors for the following code. I really appreciate your help. fit0<-lm(Y~1, data= mydata) fit.final<- lm(Y~X1+X2+X3+.....+X400, data=mydata) ??? step(fit0, scope=list(lower=fit0, upper=fit.final), data=mydata, direction="forward") step(fit.final, scope=list(lower=fit.final, upper=fit0), data=mydata, direction="backward")   Best,Farnoosh S...

...is EMail message. John > #Create simple survival object > GVHDdata<-list(Time=GVHD$Time,Time30=(GVHD$Time)/30, + Age=GVHD$Age,Drug=GVHD$Drug,Died=GVHD$Died, + AgeGrp=cut(GVHD$Age,breaks=c(0,15,25,45))) > > summary(GVHD$Drug) MTX MXT+CSP 32 32 > > > > fit0<-coxph(Surv(Time30,Died)~Drug,data=GVHDdata) > summary(fit0) Call: coxph(formula = Surv(Time30, Died) ~ Drug, data = GVHDdata) n= 64 coef exp(coef) se(coef) z p Drug[T.MXT+CSP] -1.15 0.316 0.518 -2.23 0.026 exp(coef) exp(-coef) lower .95 up...

...please obey the mailing list rules and turn of First, as others have said please obey the mailing list rules and turn off html, not everyone uses an html email client. Here is your code, formatted and with line numbers added. I also fixed one error: "y" should be "status". 1. fit0 <- coxph(Surv(futime, status) ~ x1 + x2 + x3, data = data0) 2. p <- log(predict(fit0, newdata = data1, type = "expected")) 3. lp <- predict(fit0, newdata = data1, type = "lp") 4. logbase <- p - lp 5. fit1 <- glm(status ~ offset(p), family = poisson, data = data1)...

...on. The plot, which should give two lines (one for each treatment) gives only one line and a warning message. I would appreciate help getting two lines, and an explanation of the warning message. My problem may the that I have very few events in one of my strata, but I don't know. Thanks, John fit0<-coxph(Surv(Time30,Died)~strata(Rx)+Age,data=GVHDdata) plot(survfit(fit0),fun="cloglog") WARNING: Warning in xy.coords (x, y, xlabel, ylabel, log) : 2 x values <=0 omitted from logarithmic plot > print(survfit(fit0),fun="cloglog") Call: survfit.coxph(object = fit...

Here is an example, modified from the help page to use test="Cp": -------------------------------------------------------------------------------- > fit0 <- lm(sr ~ 1, data = LifeCycleSavings) > fit1 <- update(fit0, . ~ . + pop15) > fit2 <- update(fit1, . ~ . + pop75) > anova(fit0, fit1, fit2, test="Cp") Error in `[.data.frame`(table, , "Resid. Dev") : undefined columns selected > sessionInfo() R version 2....

...ta.frame") mydata[,"StepType"] <- rep(c("First","Second"),25) mydata ######################## # END Create Dataframe # ######################## ############################ # Define function to be run# ############################ DoReg <- function(x){ fit0<-lm(as.numeric(HipFlex) ~ jSex,data=x) print(summary(fit0)) cat("\nMale\n") print(contrast(fit0, list(jSex="Male"))) cat("\nFemale\n") print(contrast(fit0, list(jSex="Female"))) cat("\nDifferenc...

...- heart$start fit <- glm(transplant ~ age + surgery + year + follow, data=heart, family = binomial) heart$iptw <- ifelse(heart$transplant == 0, 1 - predict(fit, type = "response"), predict(fit, type = "response")) summary(heart$iptw) # no weights (basic calculation) fit0 <- coxph(Surv(start,stop,event)~transplant, data=heart) fit0 # fit with coxph without case-weights fit1 <- coxreg(Surv(start,stop,event)~transplant, data=heart) fit1 # fit with coxreg from eha without case-weights # coxph fit2 <- coxph(Surv(start,stop,event)~transplant + cluster(id), da...

...would I call them?? I can't use the word "response", because that's already used with a different meaning. Might "observations" be the appropriate term? ????? What do you think? ????? Thanks, ????? Spencer Graves > x0 <- c(-1, 1) > var(x0) [1] 2 > fit0 <- lm(x0~1) > vcov(fit0) ??????????? (Intercept) (Intercept)?????????? 1 > sim0 <- simulate(fit0, 10000, 1) > var(unlist(sim0)) [1] 2.006408 > x1 <- 1 > fit1 <- glm(x1~1, poisson) > coef(fit1) ?(Intercept) 4.676016e-11 > exp(coef(fit1)) (Intercept) ???????...

...2468 12.506960 4.572391 > > > ## Using the initial parameters above, > ## fit the data with a logistic curve. > para0.st <- c( + alpha=2.212, + beta=12.507/4.572, # beta in our model is xmid/scale + gamma=1/4.572 # gamma (or r) is 1/scale + ) > > fit0 <- nls( + severity~alpha/(1+exp(beta-gamma*time)), + data1, + start=para0.st, + trace=T + ) 0.1621433 : 2.2120000 2.7355643 0.2187227 0.1621427 : 2.2124095 2.7352979 0.2187056 > > ## Plot to see how the model fits the data; plot the > ## logistic cur...

Folks: I'm having trouble doing a forward variable selection using step() First, I fit an initial model: fit0 <- glm ( est~1 , data=all, subset=c(n>=25) ) then I invoke step(): fit1 <- step( fit0 , scope=list(upper=est~ pcped + pchosp + pfarm ,lower=est~1)) I get the error message: Error in eval(expr, envir, enclos) : invalid 'envir' argument I looked at the documention on ste...

Hello, I'm quite new to R and want to make a Weibull-regression with the survival package. I know how to build my "Surv"-object and how to make a standard-weibull regression with "survreg". However, I want to fit a translated or 3-parametric weibull dist to account for a failure-free time. I think I would need a new object in survreg.distributions, but I don't know how

...r <- 10 chisq.sim <- rep(NA, iter) set.seed(1) for(i in 1:iter){ s.i <- sample(subj.) # Randomize subject assignments to 'pred' groups epil3.$pred <- pred.[s.i][epil3.$subject] fit1 <- lmer(y ~ pred+(1 | subject), family = poisson, data = epil3.) fit0 <- lmer(y ~ 1+(1 | subject), family = poisson, data = epil3.) chisq.sim[i] <- anova(fit0, fit1)[2, "Chisq"] } ************simulation (MM)************ iter <- 10 chisq.sim <- rep(NA, iter) Zt <- slot(model1,"Zt") # see ?lmer-class n.grps <-...

...l, and I wonder if someone can suggest a better way to do this. The following example consists of a slight downward trend followed by a jump up after t1=4, following by a more marked downward trend after t1=5: Dat0 <- cbind(t1=1:10, x=c(1, 0, 0, 90, 99, 95, 90, 87, 80, 77)) library(mda) fit0 <- mars(Dat0[, 1, drop=FALSE], Dat0[, 2], penalty=.001) plot(Dat0, type="l") lines(Dat0[, 1], fit0$fitted.values, lty=2, col="red") Are there 'mars' options I'm missing or other software I should be using? I've got thousands...

...fferent meaning. Might "observations" be the >> appropriate term? >> >> >> ? ????? What do you think? >> ? ????? Thanks, >> ? ????? Spencer Graves >> >> >> ? > x0 <- c(-1, 1) >> ? > var(x0) >> [1] 2 >> ? > fit0 <- lm(x0~1) >> ? > vcov(fit0) >> ? ??????????? (Intercept) >> (Intercept)?????????? 1 >> ? > sim0 <- simulate(fit0, 10000, 1) >> ? > var(unlist(sim0)) >> [1] 2.006408 >> ? > x1 <- 1 >> ? > fit1 <- glm(x1~1, poisson) >>...

Dear R-helpers, I want to try to extract residuals from a multi-level linear mixed effects model, to correlate with another variable. I need to know which residuals relate to which experimental units in the lme. I can show the labels that relate to the experimental units via the command ranef(fit0)$resid which gives: 604/1/0 -1.276971e-05 604/1/1 -1.078644e-03 606/1/0 -7.391706e-03 606/1/1 8.371521e-03 610/1/0 -6.361077e-03 610/1/1 -1.090040e-03 646/1/0 -1.696881e-03 646/1/1 -6.396153e-03 But, I cannot figure out how to access the labels for each row in this table. names(unlist(ranef(fit0...

How can I extract K (number of parameters) from an AIC calculation, both to report K itself and to calculate AICc? I'm aware of the conversion from AIC -> AICc, where AICc = AIC + 2K(K+1)/(n-K-1), but not sure of how K is calculated or how to extract that value from either an AIC or logLik calculation. This is probably more of a basic statistics question than an R question, but I thank