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