Dear Michael (and all :)) Thank you very much.
I fixed my problem, I think ;) Best, RO Atenciosamente, Rosa Oliveira -- ____________________________________________________________________________ Rosa Celeste dos Santos Oliveira, E-mail: rosit...@gmail.com Tlm: +351 939355143 Linkedin: https://pt.linkedin.com/in/rosacsoliveira ____________________________________________________________________________ "Many admire, few know" Hippocrates > On 24 Sep 2015, at 08:51, Michael Dewey <li...@dewey.myzen.co.uk> wrote: > > Dear Rosa > > Please keep the list on the recipients as others may be able to help. > > See inline > > On 23/09/2015 19:19, Rosa Oliveira wrote: >> Dear Michael, >> >> *New cleaned code :) (I think :))* >> >> casedata <-read.spss("tas_05112008.sav") >> tas.data<-data.frame(casedata) >> >> #Delete patients that were not discharged >> tas.data <- tas.data[ tas.data$hosp!="si ",] >> tas.data$resultado.hosp <- ifelse(tas.data$hosp=="l", 0, 1) >> >> tas.data$tas_d2 <- >> log(ifelse(tas.data$tas_d2==88888|tas.data$tas_d2==99999, NA, >> tas.data$tas_d2)) >> tas.data$tas_d3 <- >> log(ifelse(tas.data$tas_d3==88888|tas.data$tas_d3==99999, NA, >> tas.data$tas_d3)) >> tas.data$tas_d4 <- >> log(ifelse(tas.data$tas_d4==88888|tas.data$tas_d4==99999, NA, >> tas.data$tas_d4)) >> tas.data$tas_d5 <- >> log(ifelse(tas.data$tas_d5==88888|tas.data$tas_d5==99999, NA, >> tas.data$tas_d5)) >> tas.data$tas_d6 <- >> log(ifelse(tas.data$tas_d6==88888|tas.data$tas_d6==99999, NA, >> tas.data$tas_d6)) >> >> tas.data$age <- >> ifelse(tas.data$age==88888|tas.data$age==99999, NA, tas.data$age) >> >> >> tas.data <- data.frame(tas1 = >> tas.data$tas_d2, tas2 = tas.data$tas_d3, >> tas3 = tas.data$tas_d4, >> tas4 = tas.data$tas_d5, >> tas5 = tas.data$tas_d6, >> age = tas.data$age, >> discharge = >> tas.data$resultado.hosp, id.pat=tas.data$id) >> >> # tas.data$discharge <- factor( tas.data$discharge , >> levels=c(0,1), labels = c("dead", "alive")) >> >> #select only cases that have more than 3 tas >> tas.data <- tas.data[apply(tas.data[,-8:-6], >> 1, function(x) sum(!is.na(x)))>2,] >> >> >> nsample <- n.obs <- dim(tas.data)[1] #nr of patients >> with more than 2 tas measurements >> >> tas.data.long <- data.frame( >> tas=c(unlist(tas.data[,-8:-6])), time=rep(c(0:4), each=n.obs), >> age=rep(tas.data$age, 5), discharge=rep(tas.data$discharge, 5), >> id=rep(c(1:n.obs), 5)) >> tas.data.long <- tas.data.long >> [order(tas.data.long$id),] >> >> age=tas.data$age >> >> >> >> library(verification) > > What does that do? > >> prevOR1 <- >> summary(glm(tas.data[,7]~tas.data[,4],family=binomial(link=probit))) >> ORmodel1 <- exp(prevOR1$coeff[,1])#####computes OR? >> ORmodel1 >> >> prevOR2 <- >> summary(glm(tas.data[,7]~tas.data[,4]+tas.data[,6],family=binomial(link=probit))) >> ORmodel2 <- exp(prevOR2$coeff[,1])#####computes OR? >> ORmodel2 >> > > So are you happy that those are odds ratios but you need the confidence > intervals now? > > ?confint > > may help you >> >> Nonetheless I can’t get OR confidence intervals :( and i’m not sure if I >> have it right :( >> >> Best, >> RO >> >> >> >> Atenciosamente, >> Rosa Oliveira >> >> -- >> ____________________________________________________________________________ >> >> >> Rosa Celeste dos Santos Oliveira, >> >> E-mail:rosit...@gmail.com <mailto:rosit...@gmail.com> >> Tlm: +351 939355143 >> Linkedin: https://pt.linkedin.com/in/rosacsoliveira >> ____________________________________________________________________________ >> "Many admire, few know" >> Hippocrates >> >>> On 23 Sep 2015, at 16:29, Michael Dewey <li...@dewey.myzen.co.uk >>> <mailto:li...@dewey.myzen.co.uk>> wrote: >>> >>> Dear Rosa >>> >>> Can you remove all the code which is not relevant to calculating the >>> odds ratio so we can see what is going on? >>> >>> On 23/09/2015 16:06, Rosa Oliveira wrote: >>>> Dear Michael, >>>> >>>> >>>> I found some of the errors, but others I wasn’t able to. >>>> >>>> And my huge huge problem concerns OR and OR confidence interval :( >>>> >>>> >>>> *New Corrected code:* >>>> >>>> >>>> casedata <-read.spss("tas_05112008.sav") >>>> tas.data<-data.frame(casedata) >>>> >>>> #Delete patients that were not discharged >>>> tas.data <- tas.data[ tas.data$hosp!="si ",] >>>> tas.data$resultado.hosp <- ifelse(tas.data$hosp=="l", 0, 1) >>>> >>>> tas.data$tas_d2 <- >>>> log(ifelse(tas.data$tas_d2==88888|tas.data$tas_d2==99999, NA, >>>> tas.data$tas_d2)) >>>> tas.data$tas_d3 <- >>>> log(ifelse(tas.data$tas_d3==88888|tas.data$tas_d3==99999, NA, >>>> tas.data$tas_d3)) >>>> tas.data$tas_d4 <- >>>> log(ifelse(tas.data$tas_d4==88888|tas.data$tas_d4==99999, NA, >>>> tas.data$tas_d4)) >>>> tas.data$tas_d5 <- >>>> log(ifelse(tas.data$tas_d5==88888|tas.data$tas_d5==99999, NA, >>>> tas.data$tas_d5)) >>>> tas.data$tas_d6 <- >>>> log(ifelse(tas.data$tas_d6==88888|tas.data$tas_d6==99999, NA, >>>> tas.data$tas_d6)) >>>> >>>> tas.data$age <- >>>> ifelse(tas.data$age==88888|tas.data$age==99999, NA, tas.data$age) >>>> >>>> >>>> tas.data <- data.frame(tas1 = >>>> tas.data$tas_d2, tas2 = tas.data$tas_d3, >>>> tas3 = tas.data$tas_d4, >>>> tas4 = tas.data$tas_d5, >>>> tas5 = tas.data$tas_d6, >>>> age = tas.data$age, >>>> discharge = >>>> tas.data$resultado.hosp, id.pat=tas.data$id) >>>> >>>> # tas.data$discharge <- factor( tas.data$discharge , >>>> levels=c(0,1), labels = c("dead", "alive")) >>>> >>>> #select only cases that have more than 3 tas >>>> tas.data <- tas.data[apply(tas.data[,-8:-6], >>>> 1, function(x) sum(!is.na(x)))>2,] >>>> >>>> >>>> nsample <- n.obs <- dim(tas.data)[1] #nr of patients >>>> with more than 2 tas measurements >>>> >>>> tas.data.long <- data.frame( >>>> tas=c(unlist(tas.data[,-8:-6])), time=rep(c(0:4), each=n.obs), >>>> age=rep(tas.data$age, 5), discharge=rep(tas.data$discharge, 5), >>>> id=rep(c(1:n.obs), 5)) >>>> tas.data.long <- tas.data.long >>>> [order(tas.data.long$id),] >>>> >>>> age=tas.data$age >>>> >>>> ################################################################################################## >>>> #PLOT EMPIRICAL MEANS OF CRP FOR ALIVE DEATh >>>> ################################################################################################## >>>> mean.alive <- >>>> apply(tas.data[tas.data$discharge==0, -8:-6], 2, mean, na.rm=T) >>>> mean.dead <- >>>> apply(tas.data[tas.data$discharge==1, -8:-6], 2, mean, na.rm=T) >>>> stderr <- function(x) >>>> sqrt(var(x,na.rm=TRUE)/length(na.omit(x))) >>>> se.alive <- >>>> apply(tas.data[tas.data$discharge==0, -8:-6], 2, stderr) >>>> se.dead <- >>>> apply(tas.data[tas.data$discharge==1, -8:-6], 2, stderr) >>>> summarytas <- data.frame (media = c(mean.dead, >>>> mean.alive), >>>> standarderror = c( se.dead, >>>> se.alive), discharge = c(rep("dead",5), rep("alive", 5))) >>>> >>>> >>>> ggplot(summarytas, aes(x=c(c(1:5), c(1:5)), y=mean, colour=discharge)) + >>>> geom_errorbar(aes(ymin=mean-2* standarderror, ymax=media+2* >>>> standarderror), width=.1) + >>>> scale_color_manual(values=c("blue", "red")) + >>>> theme(legend.text=element_text(size=20), >>>> axis.text=element_text(size=16), >>>> axis.title=element_text(size=20,face="bold")) + >>>> scale_x_continuous(name="Days") + >>>> scale_y_continuous(name="log tas") + >>>> geom_line() + >>>> geom_point() >>>> >>>> >>>> library(verification) >>>> prev <- >>>> summary(glm(tas.data[,7]~tas.data[,4],family=binomial(link=probit))) >>>> answer = c(prev$coefficients[,1:2]) >>>> >>>> >>>> roc.plot(tas.data[,7], prev, show.thres = FALSE, legen=F ) >>>> >>>> >>>> >>>> modelo<-function (dataainit) >>>> { >>>> >>>> #dataa<-tas.data >>>> dataa<-dataainit >>>> >>>> dataa$ident<-seq(1:90) >>>> tas.6days<-cbind(id=rep(dataa$id,5),tas=c(dataa$tas_d2, dataa$tas_d3, >>>> dataa$tas_d4, dataa$tas_d5, dataa$tas_d6), >>>> time=rep(c(2:6)-2, each=90), out.come=rep(dataa$hosp,5), >>>> ident=rep(dataa$ident,5)) >>>> tas.6days<-data.frame(tas.6days[order(tas.6days[,1]),]) >>>> tas.6days$tas[tas.6days$tas==88888|tas.6days$tas==99999 ]<-NA >>>> #mixed model for the longitudinal tas >>>> lme.1 <- lme(tas~ time+1, random = ~ time+1 |ident, data=tas.6days, >>>> na.action=na.exclude ) >>>> #Random intercept and slopes >>>> pred.lme<-predict(lme.1) >>>> lme.intercept<-lme.1$coef$random$ident[,1]+lme.1$coef$fixed[1] >>>> lme.slope<- lme.1$coef$random$ident[,2]+lme.1$coef$fixed[2] >>>> selector<-as.numeric(names(lme.intercept)) #to select not NA rows. Apply >>>> to the vector in the dataset >>>> test(dataa$intercept[resultado.hosp==1], >>>> dataa$intercept[resultado.hosp==0]) >>>> print(summary(model.surv1)) >>>> return(model.surv1$coef) >>>> } >>>> >>>> >>>> *CONSOLE RESULT: (errors in red)* >>>> >>>> > casedata <-read.spss("tas_05112008.sav") >>>> Warning message: >>>> In read.spss("tas_05112008.sav") : >>>> tas_05112008.sav: Unrecognized record type 7, subtype 18 encountered >>>> in system file >>>> > tas.data<-data.frame(casedata) >>>> > >>>> > #Delete patients that were not discharged >>>> > tas.data <- tas.data[ tas.data$hosp!="si ",] >>>> > tas.data$resultado.hosp <- ifelse(tas.data$hosp=="l", 0, 1) >>>> > >>>> > tas.data$tas_d2 <- >>>> log(ifelse(tas.data$tas_d2==88888|tas.data$tas_d2==99999, NA, >>>> tas.data$tas_d2)) >>>> > tas.data$tas_d3 <- >>>> log(ifelse(tas.data$tas_d3==88888|tas.data$tas_d3==99999, NA, >>>> tas.data$tas_d3)) >>>> > tas.data$tas_d4 <- >>>> log(ifelse(tas.data$tas_d4==88888|tas.data$tas_d4==99999, NA, >>>> tas.data$tas_d4)) >>>> > tas.data$tas_d5 <- >>>> log(ifelse(tas.data$tas_d5==88888|tas.data$tas_d5==99999, NA, >>>> tas.data$tas_d5)) >>>> > tas.data$tas_d6 <- >>>> log(ifelse(tas.data$tas_d6==88888|tas.data$tas_d6==99999, NA, >>>> tas.data$tas_d6)) >>>> > >>>> > tas.data$age <- >>>> ifelse(tas.data$age==88888|tas.data$age==99999, NA, tas.data$age) >>>> > >>>> > >>>> > tas.data <- data.frame(tas1 = >>>> tas.data$tas_d2, tas2 = tas.data$tas_d3, >>>> + tas3 = tas.data$tas_d4, >>>> tas4 = tas.data$tas_d5, >>>> + tas5 = tas.data$tas_d6, >>>> age = tas.data$age, >>>> + discharge = >>>> tas.data$resultado.hosp, id.pat=tas.data$id) >>>> > >>>> > # tas.data$discharge <- factor( tas.data$discharge >>>> , levels=c(0,1), labels = c("dead", "alive")) >>>> > >>>> > #select only cases that have more than 3 tas >>>> > tas.data <- tas.data[apply(tas.data[,-8:-6], >>>> 1, function(x) sum(!is.na(x)))>2,] >>>> > >>>> > >>>> > >>>> > nsample <- n.obs <- dim(tas.data)[1] #nr of >>>> patients with more than 2 tas measurements >>>> > >>>> > tas.data.long <- data.frame( >>>> tas=c(unlist(tas.data[,-8:-6])), time=rep(c(0:4), each=n.obs), >>>> age=rep(tas.data$age, 5), discharge=rep(tas.data$discharge, 5), >>>> + id=rep(c(1:n.obs), 5)) >>>> > tas.data.long <- tas.data.long >>>> [order(tas.data.long$id),] >>>> > >>>> > age=tas.data$age >>>> > >>>> > >>>> ################################################################################################## >>>> > #PLOT EMPIRICAL MEANS OF CRP FOR ALIVE DEATh >>>> > >>>> ################################################################################################## >>>> > mean.alive <- >>>> apply(tas.data[tas.data$discharge==0, -8:-6], 2, mean, na.rm=T) >>>> > mean.dead <- >>>> apply(tas.data[tas.data$discharge==1, -8:-6], 2, mean, na.rm=T) >>>> > stderr <- function(x) >>>> sqrt(var(x,na.rm=TRUE)/length(na.omit(x))) >>>> > se.alive <- >>>> apply(tas.data[tas.data$discharge==0, -8:-6], 2, stderr) >>>> > se.dead <- >>>> apply(tas.data[tas.data$discharge==1, -8:-6], 2, stderr) >>>> > summarytas <- data.frame (media = c(mean.dead, >>>> mean.alive), >>>> + standarderror = c( se.dead, >>>> se.alive), discharge = c(rep("dead",5), rep("alive", 5))) >>>> > >>>> > >>>> > ggplot(summarytas, aes(x=c(c(1:5), c(1:5)), y=mean, >>>> colour=discharge)) + >>>> + geom_errorbar(aes(ymin=mean-2* standarderror, ymax=media+2* >>>> standarderror), width=.1) + >>>> + scale_color_manual(values=c("blue", "red")) + >>>> + theme(legend.text=element_text(size=20), >>>> axis.text=element_text(size=16), >>>> axis.title=element_text(size=20,face="bold")) + >>>> + scale_x_continuous(name="Days") + >>>> + scale_y_continuous(name="log tas") + >>>> + geom_line() + >>>> + geom_point() >>>> Error in mean - 2 * standarderror : >>>> non-numeric argument to binary operator >>>> > >>>> > >>>> > library(verification) >>>> > prev <- >>>> summary(glm(tas.data[,7]~tas.data[,4],family=binomial(link=probit))) >>>> > answer = c(prev$coefficients[,1:2]) >>>> > >>>> > >>>> > roc.plot(tas.data[,7], prev, show.thres = FALSE, legen=F ) >>>> Error in is.finite(x) : default method not implemented for type 'list' >>>> > >>>> > >>>> > >>>> > modelo<-function (dataainit) >>>> + >>>> + { >>>> + >>>> + #dataa<-tas.data >>>> + dataa<-dataainit >>>> + >>>> + dataa$ident<-seq(1:90) >>>> + tas.6days<-cbind(id=rep(dataa$id,5),tas=c(dataa$tas_d2, dataa$tas_d3, >>>> + dataa$tas_d4, dataa$tas_d5, dataa$tas_d6), >>>> + time=rep(c(2:6)-2, each=90), out.come=rep(dataa$hosp,5), >>>> ident=rep(dataa$ident,5)) >>>> + >>>> + tas.6days<-data.frame(tas.6days[order(tas.6days[,1]),]) >>>> + tas.6days$tas[tas.6days$tas==88888|tas.6days$tas==99999 ]<-NA >>>> + >>>> + #mixed model for the longitudinal tas >>>> + lme.1 <- lme(tas~ time+1, random = ~ time+1 |ident, data=tas.6days, >>>> na.action=na.exclude ) >>>> + >>>> + #Random intercept and slopes >>>> + pred.lme<-predict(lme.1) >>>> + lme.intercept<-lme.1$coef$random$ident[,1]+lme.1$coef$fixed[1] >>>> + lme.slope<- lme.1$coef$random$ident[,2]+lme.1$coef$fixed[2] >>>> + selector<-as.numeric(names(lme.intercept)) #to select not NA rows. >>>> Apply to the vector in the dataset >>>> + >>>> + test(dataa$intercept[resultado.hosp==1], >>>> dataa$intercept[resultado.hosp==0]) >>>> + >>>> + print(summary(model.surv1)) >>>> + return(model.surv1$coef) >>>> + >>>> + } >>>> > >>>> >>>> I can’t get the OR and OR CI :( >>>> >>>> >>>> *DATA:* >>>> >>>> >>>> >>>> >>>> >>>> >>>> Best, >>>> >>>> RO >>>> >>>> >>>> >>>> >>>> Atenciosamente, >>>> Rosa Oliveira >>>> >>>> -- >>>> ____________________________________________________________________________ >>>> >>>> >>>> >>>> >>>> Rosa Celeste dos Santos Oliveira, >>>> >>>> E-mail:rosit...@gmail.com <http://gmail.com> <mailto:rosit...@gmail.com> >>>> Tlm: +351 939355143 >>>> Linkedin: https://pt.linkedin.com/in/rosacsoliveira >>>> ____________________________________________________________________________ >>>> "Many admire, few know" >>>> Hippocrates >>>> >>>>> On 23 Sep 2015, at 12:02, Michael Dewey <li...@dewey.myzen.co.uk >>>>> <mailto:li...@dewey.myzen.co.uk> >>>>> <mailto:li...@dewey.myzen.co.uk>> wrote: >>>>> >>>>> Dear Rosa >>>>> >>>>> It would help if you posted the error messages where they occur so >>>>> that we can see which of your commands caused which error. However see >>>>> comment inline below. >>>>> >>>>> On 22/09/2015 22:17, Rosa Oliveira wrote: >>>>>> Dear all, >>>>>> >>>>>> >>>>>> I’m trying to compute Odds ratio and OR confidence interval. >>>>>> >>>>>> I’m really naive, sorry for that. >>>>>> >>>>>> >>>>>> I attach my data and my code. >>>>>> >>>>>> I’m having lots of errors: >>>>>> >>>>>> 1. Error in data.frame(tas1 = tas.data$tas_d2, tas2 = >>>>>> tas.data$tas_d3, tas3 = tas.data$tas_d4, : >>>>>> arguments imply differing number of rows: 90, 0 >>>>>> >>>>> >>>>> At least one of tas_d2, tas_d3, tas_d4 does not exist >>>>> >>>>> I suggest fixing that one and hoping the rest go away. >>>>> >>>>>> 2. Error in data.frame(tas = c(unlist(tas.data[, -8:-6])), time = >>>>>> rep(c(0:4), : >>>>>> arguments imply differing number of rows: 630, 450, 0 >>>>>> >>>>>> 3. Error: object 'tas.data.long' not found >>>>>> >>>>>> 4. Error in data.frame(media = c(mean.dead, mean.alive), >>>>>> standarderror = c(se.dead, : >>>>>> arguments imply differing number of rows: 14, 10 >>>>>> >>>>>> 5. Error in ggplot(summarytas, aes(x = c(c(1:5), c(1:5)), y = mean, >>>>>> colour = discharge)) : >>>>>> object 'summarytas' not found >>>>>> >>>>>> 6. Error in summary(glm(tas.data[, 6] ~ tas.data[, 4], family = >>>>>> binomial(link = probit))) : >>>>>> error in evaluating the argument 'object' in selecting a method for >>>>>> function 'summary': Error in eval(expr, envir, enclos) : y values >>>>>> must be 0 <= y <= 1 >>>>>> >>>>>> 7. Error in wilcox.test.default(pred[obs == 1], pred[obs == 0], >>>>>> alternative = "great") : >>>>>> not enough (finite) 'x' observations >>>>>> In addition: Warning message: >>>>>> In is.finite(x) & apply(pred, 1, f) : >>>>>> longer object length is not a multiple of shorter object length >>>>>> >>>>>> >>>>>> and off course I’m not getting OR. >>>>>> >>>>>> Nonetheless all this errors, I think I have not been able to compute >>>>>> de code to get OR and OR confidence interval. >>>>>> >>>>>> >>>>>> Can anyone help me please. It’s really urgent. >>>>>> >>>>>> PLEASE >>>>>> >>>>>> THE CODE: >>>>>> >>>>>> the hospital outcome is discharge. >>>>>> >>>>>> require(gdata) >>>>>> library(foreign) >>>>>> library(nlme) >>>>>> library(lme4) >>>>>> library(boot) >>>>>> library(MASS) >>>>>> library(Hmisc) >>>>>> library(plotrix) >>>>>> library(verification) >>>>>> library(mvtnorm) >>>>>> library(statmod) >>>>>> library(epiR) >>>>>> >>>>>> ######################################################################################### >>>>>> # Data preparation >>>>>> # >>>>>> ######################################################################################### >>>>>> >>>>>> setwd("/Users/RO/Desktop") >>>>>> >>>>>> casedata <-read.spss("tas_05112008.sav") >>>>>> tas.data<-data.frame(casedata) >>>>>> >>>>>> #Delete patients that were not discharged >>>>>> tas.data <- tas.data[ tas.data$hosp!="si ",] >>>>>> tas.data$resultado.hosp <- ifelse(tas.data$hosp=="l", 0, 1) >>>>>> >>>>>> tas.data$tas_d2 <- >>>>>> log(ifelse(tas.data$tas_d2==88888|tas.data$tas_d2==99999, NA, >>>>>> tas.data$tas_d2)) >>>>>> tas.data$tas_d3 <- >>>>>> log(ifelse(tas.data$tas_d3==88888|tas.data$tas_d3==99999, NA, >>>>>> tas.data$tas_d3)) >>>>>> tas.data$tas_d4 <- >>>>>> log(ifelse(tas.data$tas_d4==88888|tas.data$tas_d4==99999, NA, >>>>>> tas.data$tas_d4)) >>>>>> tas.data$tas_d5 <- >>>>>> log(ifelse(tas.data$tas_d5==88888|tas.data$tas_d5==99999, NA, >>>>>> tas.data$tas_d5)) >>>>>> tas.data$tas_d6 <- >>>>>> log(ifelse(tas.data$tas_d6==88888|tas.data$tas_d6==99999, NA, >>>>>> tas.data$tas_d6)) >>>>>> >>>>>> tas.data$age <- >>>>>> ifelse(tas.data$age==88888|tas.data$age==99999, NA, tas.data$age) >>>>>> >>>>>> >>>>>> tas.data <- data.frame(tas1 = >>>>>> tas.data$tas_d2, tas2 = tas.data$tas_d3, >>>>>> tas3 = tas.data$tas_d4, >>>>>> tas4 = tas.data$tas_d5, >>>>>> tas5 = tas.data$tas_d6, >>>>>> age = tas.data$age, >>>>>> discharge = >>>>>> tas.data$resultado.hosp, id.pat=tas.data$ID) >>>>>> >>>>>> # tas.data$discharge <- factor( tas.data$discharge >>>>>> , levels=c(0,1), labels = c("dead", "alive")) >>>>>> >>>>>> #select only cases that have more than 3 tas >>>>>> tas.data <- tas.data[apply(tas.data[,-8:-6], >>>>>> 1, function(x) sum(!is.na(x)))>2,] >>>>>> >>>>>> >>>>>> >>>>>> nsample <- n.obs <- dim(tas.data)[1] #nr of >>>>>> patients with more than 2 tas measurements >>>>>> >>>>>> tas.data.long <- data.frame( >>>>>> tas=c(unlist(tas.data[,-8:-6])), time=rep(c(0:4), each=n.obs), >>>>>> age=rep(tas.data$age, 5), discharge=rep(tas.data$discharge, 5), >>>>>> id=rep(c(1:n.obs), 5)) >>>>>> tas.data.long <- tas.data.long >>>>>> [order(tas.data.long$id),] >>>>>> >>>>>> age=tas.data$age >>>>>> >>>>>> ################################################################################################## >>>>>> #PLOT EMPIRICAL MEANS OF CRP FOR ALIVE DEATh >>>>>> ################################################################################################## >>>>>> mean.alive <- >>>>>> apply(tas.data[tas.data$discharge==0, -8:-6], 2, mean, na.rm=T) >>>>>> mean.dead <- >>>>>> apply(tas.data[tas.data$discharge==1, -8:-6], 2, mean, na.rm=T) >>>>>> stderr <- function(x) >>>>>> sqrt(var(x,na.rm=TRUE)/length(na.omit(x))) >>>>>> se.alive <- >>>>>> apply(tas.data[tas.data$discharge==0, -8:-6], 2, stderr) >>>>>> se.dead <- >>>>>> apply(tas.data[tas.data$discharge==1, -8:-6], 2, stderr) >>>>>> summarytas <- data.frame (media = c(mean.dead, >>>>>> mean.alive), >>>>>> standarderror = c( se.dead, >>>>>> se.alive), discharge = c(rep("dead",5), rep("alive", 5))) >>>>>> >>>>>> >>>>>> ggplot(summarytas, aes(x=c(c(1:5), c(1:5)), y=mean, >>>>>> colour=discharge)) + >>>>>> geom_errorbar(aes(ymin=mean-2* standarderror, ymax=media+2* >>>>>> standarderror), width=.1) + >>>>>> scale_color_manual(values=c("blue", "red")) + >>>>>> theme(legend.text=element_text(size=20), >>>>>> axis.text=element_text(size=16), >>>>>> axis.title=element_text(size=20,face="bold")) + >>>>>> scale_x_continuous(name="Days") + >>>>>> scale_y_continuous(name="log tas") + >>>>>> geom_line() + >>>>>> geom_point() >>>>>> >>>>>> >>>>>> library(verification) >>>>>> prev <- >>>>>> summary(glm(tas.data[,6]~tas.data[,4],family=binomial(link=probit))) >>>>>> answer = c(prev$coefficients[,1:2]) >>>>>> >>>>>> >>>>>> roc.plot(tas.data[,6], prev, show.thres = FALSE, legen=F ) >>>>>> >>>>>> >>>>>> modelo<-function (dataainit) >>>>>> >>>>>> { >>>>>> >>>>>> #dataa<-tas.data >>>>>> dataa<-dataainit >>>>>> >>>>>> dataa$ident<-seq(1:90) >>>>>> tas.6days<-cbind(id=rep(dataa$id,5),tas=c(dataa$tas_d2, dataa$tas_d3, >>>>>> dataa$tas_d4, dataa$tas_d5, dataa$tas_d6), >>>>>> time=rep(c(2:6)-2, each=90), out.come=rep(dataa$hosp,5), >>>>>> ident=rep(dataa$ident,5)) >>>>>> >>>>>> tas.6days<-data.frame(tas.6days[order(tas.6days[,1]),]) >>>>>> tas.6days$tas[tas.6days$tas==88888|tas.6days$tas==99999 ]<-NA >>>>>> >>>>>> #mixed model for the longitudinal tas >>>>>> lme.1 <- lme(tas~ time+1, random = ~ time+1 |ident, data=tas.6days, >>>>>> na.action=na.exclude ) >>>>>> >>>>>> #Random intercept and slopes >>>>>> pred.lme<-predict(lme.1) >>>>>> lme.intercept<-lme.1$coef$random$ident[,1]+lme.1$coef$fixed[1] >>>>>> lme.slope<- lme.1$coef$random$ident[,2]+lme.1$coef$fixed[2] >>>>>> selector<-as.numeric(names(lme.intercept)) #to select not NA rows. >>>>>> Apply to the vector in the dataset >>>>>> >>>>>> test(dataa$intercept[resultado.hosp==1], >>>>>> dataa$intercept[resultado.hosp==0]) >>>>>> >>>>>> print(summary(model.surv1)) >>>>>> return(model.surv1$coef) >>>>>> >>>>>> } >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> Best, >>>>>> RO >>>>>> >>>>>> >>>>>> >>>>>> ______________________________________________ >>>>>> R-help@r-project.org <mailto:R-help@r-project.org> >>>>>> <mailto:R-help@r-project.org> mailing list -- To >>>>>> UNSUBSCRIBE and more, see >>>>>> 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. >>>>>> >>>>> >>>>> -- >>>>> Michael >>>>> http://www.dewey.myzen.co.uk/home.html >>>> >>> >>> -- >>> Michael >>> http://www.dewey.myzen.co.uk/home.html >> > > -- > Michael > http://www.dewey.myzen.co.uk/home.html ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
