library(survival)
attach(ghn) 
     # puts the variables into workspace, so you can use 'X6.MP.time' directly
dead <- c(X6.MP.dead, control.dead)
time <- c(X6.MP.time, control.time) 
group <- c(rep(1,21),rep(0,21))

################### function survreg()
  ghn.wei <- survreg(Surv(time, dead) ~ group, dist="wei")
names(ghn.wei)  # this tells us the following names in the object ghn.exp: 
# "coefficients"      "icoef"             "var"              
# "loglik"            "iter"              "linear.predictors"
# "df"                "scale"             "idf"              
# "df.residual"       "terms"             "means"            
# "call"              "dist"              "y"    
summary(ghn.wei) # prints out some of the results, such as the loglik. 
# "call" is standard, it give the function used to create the object (ghn.wei)
 ghn.wei$dist #  "weibull" - says what distribution was used 

 
ghn.wei$coef # gives the coeffs. (notice we don't need the full "coefficients")


summary(ghn.wei) # gives output which includes this: 
#             Value Std. Error     z        p
#(Intercept)  2.248      0.166 13.55 8.30e-42
#group        1.267      0.311  4.08 4.51e-05
#Log(scale)  -0.312      0.147 -2.12 3.43e-02
#
# Scale= 0.732 

# the Std. Error values comes from the square roots of the variance-covariance matrix.
ghn.wei$var # gives the variance-covariance matrix. 
 #           (Intercept)       group  Log(scale)
#(Intercept)  0.02754664 -0.03032506 -0.00661610
#group       -0.03032506  0.09649717  0.01572686
#Log(scale)  -0.00661610  0.01572686  0.02169490

sqrt( 0.02754664) # is  0.1659718, which rounds to 0.166, the value given by summary()
 
ghn.wei$scale #  is the scale, exp(-0.312)= 0.732  (denoted sigma or 1/alpha, in my notes) 

ghnwei$loglik # gives loglikelihood for the fitted model, and the model with only the intercept
ghnwei$linear.pred # gives the result intercept + 1.267 times group = 1
                   # you probable have to look at the values to follow this. 
ghn.wei$df  # the number of degrees of freedom in the model. 
ghn.wei$df.residual# the number of degrees of freedom of the residual (not discussed yet.)

########## to get fitted values, use function predict(), e.g. 

 medianW  <- predict(ghn.wei,newdata=list(group=0:1),type="uquantile",p=0.5,se.fit=T)
   # newdata=list(group=0:1) is used to get predicted values for group=0 and group=1
  names(medianW) # shows there are two parts to the object: 
 #  "fit"    "se.fit"   and [1] 
exp(medianW$fit[1]) # gives the median for groups 0.  7.24
exp(medianW$fit[2]) # gives the median for groups 1.  25.72
 # to get a 95%  conf. int, use
round(cbind( exp(medianW$fit[1]-1.96*medianW$se.fit[1]), exp(medianW$fit[1]+1.96*medianW$se.fit[1]) ),2 )
 #  this gives: 5.08 10.32. round(x,y) rounds x to y decimal places. 
 lquartW  <- predict(ghn.wei,newdata=list(group=0:1),type="uquantile",p=0.25,se.fit=T)
  #  this gives the lower quartile. How would you get the upper quartile? 

########## to find estimate probability of surviving to 10 weeks:
  'arg' should be one of “response”, “link”, “lp”, “linear”, “terms”, “quantile”, “uquantile”
 predict(ghn.wei,type="link") just give predictor on log scale mu + beta : 2.248 , 3.5 
link and lp are the same, as is linear 
 predict(ghn.wei,type="response") two value, 21 each: 33.639028-exp(3.5) 9.472116=exp(2.248)
terms - don't know. 
 predict(ghn.wei,type="quantile",p=0.5,se.fit=T) # on original time  scale, 
  # standard error estimate sensibel on linear predictor scale. 
 whereas 'uquantile' is on linear predictor scale. 
 
newdata=list(0:10 link j

# the survival function for the Weibull is 

 a <- predict(ghn.wei,newdata=list(0:1),type="quantile")
newdata=list(0:1), type="")
   # newdata=list(group=0:1) is used to get predicted values for group=0 and group=1
 a <- predict(ghn.wei,newdata=list(group=0:1),type="quantile",p=0.5,se.fit=T)

 predict(ghn.wei,newdata=list(group=0:1),type="quantile",p=0.4)
 predict(ghn.wei,newdata=list(group=0:1),type="uquantile",p=0.5)