[R] doubt with Odds ratio - URGENT HELP NEEDED

[Rosa Oliveira]

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

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