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
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>
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 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
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 <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>> 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> 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
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
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: rosita21 at 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 <lists at 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: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
Dear Rosa, coefficents of a probit-regression do not have a odds-ratio interpretation, you should use a logit link for that. cheers. Am 24.09.2015 um 09:51 schrieb Michael Dewey:> 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 >> >-- Eik Vettorazzi Department of Medical Biometry and Epidemiology University Medical Center Hamburg-Eppendorf Martinistr. 52 20246 Hamburg T ++49/40/7410-58243 F ++49/40/7410-57790 -- _____________________________________________________________________ Universit?tsklinikum Hamburg-Eppendorf; K?rperschaft des ?ffentlichen Rechts; Gerichtsstand: Hamburg | www.uke.de Vorstandsmitglieder: Prof. Dr. Burkhard G?ke (Vorsitzender), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Pr?l?, Rainer Schoppik _____________________________________________________________________ SAVE PAPER - THINK BEFORE PRINTING