#buysetest3b
# weibull 
# updated 01/06/04 to accommodate censoring by data augmentation

model 
{

for(i in 1:N) {
	lambda[i] ~ dexp(1)
	logmu[i] <-  beta[1] + beta[2] * trt[i] + beta[3] * strat[i]
	mu[i] <-  exp( logmu[i])
		
	#z[i] <- surv[i] / (mu[i] * sqrt(2 * lambda[i]) )							# z is standard half normal
	z[i] <- pow(surv[i] / mu[i], alpha) / sqrt(2 * lambda[i]) 							# z is standard half normal
	
	zeroes[i] <- 0
	psi[i] <- -log(alpha) - (alpha - 1.0) * log(surv[i]) + log(2 * lambda[i]) /2 + alpha * logmu[i]  + pow( z[i] , 2 ) / 2 + log(2 * 3.14159265) / 2 + 2	# minus log likelihood
	zeroes[i] ~ dpois( psi[i])
	zeroes2[i] <- 0
	#mark[i] ~ dnorm( mumark[i], deltastar1 )
	mumark[i] <- betas[1] + betas[2] * trt[i] + betas[3] * strat[i] + deltastar2 * z[i]
	s[i] ~ dnorm( mumark[i], deltastar1)
	#psi2[i] <-  -log( deltastar1)/2  + log(2 * 3.14159265) /2 + deltastar1 * pow( s[i] - mumark[i], 2) /2
	#zeroes2[i] ~ dpois( psi2[i] )
	surv[i] ~ dunif(cent[i], 2000)
}	

for (i in 1:N) {
 	surv2[i] ~ dweib( alpha2, mu2[i] )I(cent2[i],)
	 mu2[i] <-  exp( -(beta2[1] + beta2[2] * trt[i] + beta2[3] * strat[i]))
}


for (j in 1:3) {
	beta[j] ~ dflat()
	beta2[j] ~ dflat()
	betas[j] ~ dflat()
	beta2star[j] <- beta2[j] / alpha2
}
	alpha ~ dexp(0.1)
	alpha2 ~ dexp(0.1)
	deltastar3 <- 1
	at <- nu / 2
	bt1 <- R[1,1] - pow(R[1,2], 2) / R[2,2]
	bt <- bt1 / 2
	deltastar1 ~ dgamma(at, bt)
	priorprec <- R[2,2] * deltastar1
	priormean <- R[1,2] / R[2,2]
	deltastar2 ~ dnorm(priormean, priorprec)
	
	Sigma[2,2] <- 1 / deltastar3
	Sigma[1,2] <- deltastar2 * Sigma[2,2]
	Sigma[1,1] <- 1 / deltastar1 + pow(Sigma[1,2],2) / Sigma[2,2]
	
	gammaz <- Sigma[1,2] / sqrt( Sigma[1,1] * Sigma[2,2])
	RE3 <- beta[2] / ( betas[2] / sqrt( Sigma[1,1]) )
}

ACTG 175 crna
data
list( R = structure(.Data = c(0.25, 0, 0, 1), .Dim = c(2,2)), nu = 4)

inits
list(beta = c(4.865, -0.444, -1.655), alpha = 0.911, alpha2 = 0.911, beta2 = 
c(5, -0.4, -1.6), deltastar1 = 4, betas = c(-0.528, -0.9647, -0.5948), deltastar2 = 0.5)
list(beta = c(5.637, 0.44, -0.587), alpha = 1.850, alpha2 = 1.850, beta2 = 
c(7, 0.7, -0.9), deltastar1 = 1, betas = c(-0.2472, -0.5951, -0.1568), deltastar2 = 0)
list(beta = c(6.409, 1.324, 0.481), alpha = 3.755, alpha2 = 3.755, beta2 = 
c(19, 4, 1), deltastar1 = 9, betas = c(0.0336, -0.2255, 0.3812), deltastar2 = -0.5)

 calcinits( 5.637, 0.193, c(-4,0,4))
[1] 4.865
[1] 5.637
[1] 6.409
> calcinits(0.440, 0.221, c(-4,0,4))
[1] -0.444
[1] 0.44
[1] 1.324
> calcinits(-0.587, 0.267, c(-4,0,4))
[1] -1.655
[1] -0.587
[1] 0.481

[1] 3.7546685 1.8496566 0.9111935

Call:
lm(formula = crna ~ trt + strat)

Residuals:
     Min       1Q   Median       3Q      Max
-1.72582 -0.39482 -0.00368  0.38232  1.54918

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.2472     0.0702  -3.521 0.000583 ***
trt          -0.5951     0.0924  -6.441 1.83e-09 ***
strat        -0.1568     0.1345  -1.166 0.245702
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Residual standard error: 0.5444 on 138 degrees of freedom
Multiple R-Squared: 0.2416,     Adjusted R-squared: 0.2306
F-statistic: 21.98 on 2 and 138 DF,  p-value: 5.185e-09

>  calcinits( -0.2472, 0.0702, c(-4,0,4))
[1] -0.528
[1] -0.2472
[1] 0.0336
> calcinits( -0.5951, 0.0924, c(-4,0,4))
[1] -0.9647
[1] -0.5951
[1] -0.2255
> calcinits( -0.1568, 0.1345, c(-4,0,4))
[1] -0.6948
[1] -0.1568
[1] 0.3812


Call:
survreg(formula = Surv(surv, fail) ~ trt + strat, dist = "weibull")
             Value Std. Error     z         p
(Intercept)  5.637      0.193 29.16 5.68e-187
trt          0.440      0.221  1.99  4.67e-02
strat       -0.587      0.267 -2.20  2.80e-02
Log(scale)  -0.615      0.177 -3.48  5.07e-04

Scale= 0.541

Weibull distribution
Loglik(model)= -194.4   Loglik(intercept only)= -198.2
        Chisq= 7.44 on 2 degrees of freedom, p= 0.024
Number of Newton-Raphson Iterations: 8
n= 141


for ACTG152 p24  use a152buysep24b.dat
data
list(N = 161, R = structure(.Data = c(0.25, 0, 0, 1), .Dim = c(2,2)), nu = 4)

inits for weibull

list(beta = c(3.43, 0.45, 0.30), alpha = 1.54, alpha2 = 1.54, beta2 = 
c(3.43, 
0.30, 0.45), betas = c(-0.409, -0.350, 0.10), deltastar1 = 5, 
deltastar2 = -0.5)



Node statistics
	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	Sigma[1,1]	0.3253	0.05453	0.001708	0.2439	0.3166	0.458	4001	40000
	Sigma[1,2]	-0.133	0.1613	0.007421	-0.4382	-0.1375	0.1904	4001	40000
	alpha	1.997	0.2349	0.006639	1.602	1.975	2.519	4001	40000
	alpha2	1.891	0.342	0.02271	1.242	1.869	2.575	4001	40000
	beta[1]	5.607	0.1432	0.007207	5.342	5.6	5.907	4001	40000
	beta[2]	0.3343	0.1673	0.008111	-0.004746	0.3412	0.6481	4001	40000
	beta[3]	-0.4826	0.2338	0.00835	-0.9228	-0.4948	-0.004066	4001	40000
	beta2star[1]	5.671	0.2099	0.009804	5.34	5.643	6.157	4001	40000
	beta2star[2]	0.4405	0.2347	0.007295	0.01497	0.4279	0.9425	4001	40000
	beta2star[3]	-0.5583	0.283	0.00688	-1.127	-0.5535	-0.008029	4001	40000
	betas[1]	-0.1404	0.1458	0.006639	-0.4384	-0.1354	0.1341	4001	40000
	betas[2]	-0.5944	0.09304	0.002002	-0.7771	-0.595	-0.4134	4001	40000
	betas[3]	-0.1576	0.1334	0.001818	-0.4161	-0.1583	0.1077	4001	40000
	gammaz	-0.224	0.2644	0.01211	-0.6647	-0.249	0.3308	4001	40000

Kernel density