[R] rms::cr.setup and Hmisc::fit.mult.impute
I have fitted a proportional odds model, but would like to compare it to a continuation ratio model. However, I am unable to fit the CR model _including_ imputated data. I guess my troubles start with settuping the data for the CR model. Any hint is appreciated! Christian library(Hmisc) library(rms) library(mice) ## simulating data (taken from rms::residuals.lrm) set.seed(1) n <- 400 age <- rnorm(n, 50, 10) blood.pressure <- rnorm(n, 120, 15) L <- .05*(age-50) + .03*(blood.pressure-120) p12 <- plogis(L) p2 <- plogis(L-1) p <- cbind(1-p12, p12-p2, p2) cp <- matrix(cumsum(t(p)) - rep(0:(n-1), rep(3,n)), byrow=TRUE, ncol=3) y <- (cp < runif(n)) %*% rep(1,3) y <- as.vector(y) ## generating missing data age[1:40]<-NA blood.pressure[30:70]<-NA ## multiple imputation using mice::mice d <- as.matrix(cbind(y, age, blood.pressure)) imp <- mice(d,seed=123) ## some cleanup rm(y, age, blood.pressure) d <-as.data.frame(d) ## proportional odds model b<-fit.mult.impute(y~age+blood.pressure, lrm, xtrans=imp, data=d) ## continuation ratio model attach(d) u <- cr.setup(y) detach(d) attach(d[u$subs,]) y <-u$y cohort<-u$cohort c <- lrm(y~cohort*(age+blood.pressure)) ## CR model with imputed data q <- fit.mult.impute(y~cohort*(age+blood.pressure), lrm, xtrans=imp) Error in model.frame.default (formula = formula, data = completed.data, : variable length differ (found for 'cohort') __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] correlated binary data and overall probability
To answer my own question: The package mvtBinaryEP was helpful! Christian Lerch schrieb: Dear All, I try to simulate correlated binary data for a clinical research project. Unfortunately, I do not come to grips with bindata(). Consider corr<-data.frame(ID=as.factor(rep(c(1:10), each=5)), task=as.factor(rep(c(1:5),10))) [this format might be more appropriate:] corr2<-data.frame(ID=as.factor(rep(c(1:10),5)), tablet=as.factor(rep(c(1:5),each=10))) Now, I want to add one column 'outcome' for the binary response: * within each 'task', the probabilty of success (1) is fixed, (say rbinom(x,1,0.7)) * within each 'ID', the outcomes are correlated (say 0.8) How can I generate this column 'outcome' with the given proporties? Many thanks for hints or even code! Regards, Christian __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] lme4 and incomplete block design
Many thanks, Bill and Emmanuel! Christian Emmanuel Charpentier schrieb: Le dimanche 08 novembre 2009 à 00:05 +0100, Christian Lerch a écrit : Dear list members, I try to simulate an incomplete block design in which every participants receives 3 out of 4 possible treatment. The outcome in binary. Assigning a binary outcome to the BIB or PBIB dataset of the package SASmixed gives the appropriate output. With the code below, fixed treatment estimates are not given for each of the 4 possible treatments, instead a kind of summary measure(?) for 'treatment' is given. block<-rep(1:24,each=3) treatment<-c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3) outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0) data<-data.frame(block,treatment,outcome) lmer(outcome~treatment +(1|block), family=binomial, data=data) Is this a problem with the recovery of interblock estimates? No... Are there special rules how incomplete block designs should look like to enable this recovery? Neither... Compare : library(lme4) Le chargement a nécessité le package : Matrix Le chargement a nécessité le package : lattice summary(lmer(outcome~treatment +(1|block), family=binomial, + data=data.frame(block<-rep(1:24,each=3), +treatment<-c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, + 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, + 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, + 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, + 4, 4, 1, 3), +outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, + 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, + 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, + 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0)) + )) Generalized linear mixed model fit by the Laplace approximation Formula: outcome ~ treatment + (1 | block) Data: data.frame(block <- rep(1:24, each = 3), treatment <- c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3), outcome <- c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0)) AIC BIC logLik deviance 86.28 93.1 -40.1480.28 Random effects: Groups NameVariance Std.Dev. block (Intercept) 0.60425 0.77734 Number of obs: 72, groups: block, 24 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.277830.71767 -1.7800.075 . treatment0.011620.25571 0.0450.964 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects: (Intr) treatment -0.892 with : summary(lmer(outcome~treatment +(1|block), family=binomial, + data=data.frame(block<-rep(1:24,each=3), +treatment<-factor(c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, +2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, +3, 2, 4, 4, 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 1, +4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, +1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3)), +outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, + 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, + 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, + 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0)) + )) Generalized linear mixed model fit by the Laplace approximation Formula: outcome ~ treatment + (1 | block) Data: data.frame(block <- rep(1:24, each = 3), treatment <- factor(c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3)), outcome <- c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1,
[R] lme4 and incomplete block design
Dear list members, I try to simulate an incomplete block design in which every participants receives 3 out of 4 possible treatment. The outcome in binary. Assigning a binary outcome to the BIB or PBIB dataset of the package SASmixed gives the appropriate output. With the code below, fixed treatment estimates are not given for each of the 4 possible treatments, instead a kind of summary measure(?) for 'treatment' is given. block<-rep(1:24,each=3) treatment<-c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3) outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0) data<-data.frame(block,treatment,outcome) lmer(outcome~treatment +(1|block), family=binomial, data=data) Is this a problem with the recovery of interblock estimates? Are there special rules how incomplete block designs should look like to enable this recovery? Any help is appreciated! Regards, Christian __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] data.frame extracting data row-wise
Nice workaround :-). Thank you. Best, Christian Original-Nachricht > Datum: Fri, 30 Oct 2009 10:54:49 +0100 > Von: Karl Ove Hufthammer > An: r-h...@stat.math.ethz.ch > Betreff: Re: [R] data.frame extracting data row-wise > On Fri, 30 Oct 2009 10:42:12 +0100 Christian Lerch > wrote: > > I am struggling with extracting data from a data frame: > > x=data.frame(a=1:11,b=100:110) > > > > What I want is a list/vector in this sence: 1 100 2 101 3 102... > > y=t(as.matrix(x)) > as.vector(y) > > Happy to help. > > -- > Karl Ove Hufthammer > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT! __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] data.frame extracting data row-wise
Dear All, I am struggling with extracting data from a data frame: x=data.frame(a=1:11,b=100:110) What I want is a list/vector in this sence: 1 100 2 101 3 102... For single rows, this works fine: as.matrix(x)[1,] For, say 2 rows, this works fine: z=c(as.matrix(x)[1,],as.matrix(x)[2,]) But z=c(as.matrix(x)[1:2,]) gives 1 2 100 101!? Is there an 'automated way' to produce a list/vector from a data frame in which data of each row are just added to the end of the previous row? Many thanks, Christian -- GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT! __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] correlated binary data and overall probability
Dear All, I try to simulate correlated binary data for a clinical research project. Unfortunately, I do not come to grips with bindata(). Consider corr<-data.frame(ID=as.factor(rep(c(1:10), each=5)), task=as.factor(rep(c(1:5),10))) [this format might be more appropriate:] corr2<-data.frame(ID=as.factor(rep(c(1:10),5)), tablet=as.factor(rep(c(1:5),each=10))) Now, I want to add one column 'outcome' for the binary response: * within each 'task', the probabilty of success (1) is fixed, (say rbinom(x,1,0.7)) * within each 'ID', the outcomes are correlated (say 0.8) How can I generate this column 'outcome' with the given proporties? Many thanks for hints or even code! Regards, Christian -- http://portal.gmx.net/de/go/dsl02 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
