regression - R survival survreg not producing a good fit -
i new using r, , trying use survival analysis in order find correlation in censored data. x data envelope mass of protostars. y data intensity of observed molecular line, , values upper limits. data is:
x <- c(17.299, 4.309, 7.368, 29.382, 1.407, 3.404, 0.450, 0.815, 1.027, 0.549, 0.018) y <- c(2.37, 0.91, 1.70, 1.97, 0.60, 1.45, 0.25, 0.16, 0.36, 0.88, 0.42) censor <- c(0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1)
i using function survreg r survival library
modeldata<-survreg(formula=surv(y,censor)~x, dist="exponential", control = list(maxiter=90))
which gives following result:
summary(modeldata) call: survreg(formula = surv(y, censor) ~ x, dist = "exponential", control = list(maxiter = 90)) value std. error z p (intercept) -0.114 0.568 -0.20 0.841 x 0.153 0.110 1.39 0.163 scale fixed @ 1 exponential distribution loglik(model)= -6.9 loglik(intercept only)= -9 chisq= 4.21 on 1 degrees of freedom, p= 0.04 number of newton-raphson iterations: 5 n= 11
however, when plot data , model using following method:
plot(x,y,pch=(censor+1)) xnew<-seq(0,30) model<-predict(modeldata,list(x=xnew)) lines(xnew,model,col="red")
i plot of x , y data; triangles censored data
i not sure going wrong. have tried different distributions, produce similar results. same true when use other data, example:
x <- c(1.14, 1.14, 1.19, 0.78, 0.43, 0.24, 0.19, 0.16, 0.17, 0.66, 0.40)
i not sure if interpreting results correctly.
i have tried other examples using same method (e.g. https://stats.idre.ucla.edu/r/examples/asa/r-applied-survival-analysis-ch-1/), , works well, far can tell.
so questions are:
am using correct function fitting data? if not, more suitable?
if correct function, why model not fitting data closely? have plotting?
thank help.
the "shape" of relationship looks concave downward, have guessed ~ log(x)
fit might be more appropriate:
dfrm <- data.frame( x = c(17.299, 4.309, 7.368, 29.382, 1.407, 3.404, 0.450, 0.815, 1.027, 0.549, 0.018), y = c(2.37, 0.91, 1.70, 1.97, 0.60, 1.45, 0.25, 0.16, 0.36, 0.88, 0.42), censor= c(0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1)) modeldata<-survreg(formula=surv(y,censor)~log(x), data=dfrm, dist="loggaussian", control = list(maxiter=90))
you code seemed appropriate:
png(); plot(y~x,pch=(censor+1),data=dfrm) xnew<-seq(0,30) model<-predict(modeldata,list(x=xnew)) lines(xnew,model,col="red"); dev.off()
modeldata call: survreg(formula = surv(y, censor) ~ log(x), data = dfrm, dist = "loggaussian", control = list(maxiter = 90)) coefficients: (intercept) log(x) 0.02092589 0.32536509 scale= 0.7861798 loglik(model)= -6.6 loglik(intercept only)= -8.8 chisq= 4.31 on 1 degrees of freedom, p= 0.038 n= 11
Comments
Post a Comment