[R] doubt with Odds ratio - URGENT HELP NEEDED
Michael Dewey
lists at dewey.myzen.co.uk
Thu Sep 24 09:51:53 CEST 2015
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:rosita21 at gmail.com <mailto:rosita21 at 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 <lists at dewey.myzen.co.uk
>> <mailto:lists at 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:rosita21 at gmail.com <http://gmail.com> <mailto:rosita21 at 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 <lists at dewey.myzen.co.uk
>>>> <mailto:lists at dewey.myzen.co.uk>
>>>> <mailto:lists at 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 at r-project.org <mailto:R-help at r-project.org>
>>>>> <mailto:R-help at 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
More information about the R-help
mailing list