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
______________________________________________
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.
