1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
-- Axiom/Fricas source file
-- Translate into .org?

-- general case
)cl p all
)pile
eta := operator 'eta
G := operator 'G
rule1 := rule (D(eta(t,beta),beta)==X; D(eta(t,beta),t)==XD*beta; D(eta(t,beta),[t,beta])==XD; eta(t,beta)==X*beta)
S := G(eta(t,beta))
H := -log(S)
h := D(H,t)
[D(expr,beta) for expr in [S,H,h]]
[rule1 expr for expr in [S,H,h]]
[rule1 D(expr,beta) for expr in [S,H,h]]

-- Proportional hazards
G(x) == exp(-exp(x))
S := G(eta(t,beta))
H := -log(S)
h := D(H,t)
D(G(x),x)
[rule1 expr for expr in [S,H,h]]
[rule1 D(expr,beta) for expr in [S,H,h]]

-- Proportional odds
G(x) == 1/(1+exp(x))
S := G(eta(t,beta))
H := -log(S)
h := D(H,t)
D(G(x),x)
[rule1 expr for expr in [S,H,h]]
[rule1 D(expr,beta) for expr in [S,H,h]]

-- Probit (is there a more canonical approach in Axiom?)
Phi := operator 'Phi
phi := operator 'phi
rule2 == rule D(Phi(x),x)==phi(x)
G(x) == Phi(-x)
S := G(eta(t,beta))
H := -log(S)
h := D(H,t)
rule2 D(G(x),x)
[rule1 rule2 expr for expr in [S,H,h]]
[rule1 rule2 D(expr,beta) for expr in [S,H,h]]

-- Additive hazards
G(x) == exp(-x)
S := G(eta(t,beta))
H := -log(S)
h := D(H,t)
D(G(x),x)
[rule1 expr for expr in [S,H,h]]
[rule1 D(expr,beta) for expr in [S,H,h]]

-- Aranda-Ordaz
G(x) == exp(-log(theta*exp(x)+1)/theta)
S := G(eta(t,beta))
H := -log(S)
h := D(H,t)
D(G(x),x)
[rule1 expr for expr in [S,H,h]]
[rule1 D(expr,beta) for expr in [S,H,h]]

-- Any other links?


-- Estimators
)cl p all

h := operator 'h
H := operator 'H
rule3 := rule (exp(logtheta)==theta; 1+theta*H(beta)==S0)

-- hazard
loghaz := log(h(beta))
rule3 D(loghaz,beta)

-- marginal hazard
marghaz := h(beta)/(1+exp(logtheta)*H(beta))
rule3 D(marghaz,beta)
rule3 D(marghaz,logtheta)
D(marghaz,beta) / (((1+exp(logtheta)*H(beta))*D(h(beta),beta)-exp(logtheta)*h(beta)*D(H(beta),beta))/(1+exp(logtheta)*H(beta))^2)
D(marghaz,logtheta) / (-exp(logtheta)*H(beta)*h(beta)/(1+exp(logtheta)*H(beta))^2)

logmarghaz := log(h(beta)/(1+exp(logtheta)*H(beta)))
rule3 D(logmarghaz,beta)
rule3 D(logmarghaz,logtheta)


-- marginal survival (identity link)
margS := (1+exp(logtheta)*H(beta))^(-exp(-logtheta)) -- identity transformation link
rule3 D(margS,beta)
rule3 D(margS,logtheta)
D(margS,beta) / (-(1+H(beta)*exp(logtheta))^(-exp(-logtheta)-1)*D(H(beta),beta))
eval(%,logtheta=1.5)
D(margS,beta) / (-margS*D(H(beta),beta)/(1+exp(logtheta)*H(beta)))
eval(%,[logtheta=1.5, H(beta)=2.0])
D(margS,logtheta) / (margS*(exp(-logtheta)*log(1+exp(logtheta)*H(beta))-H(beta)/(1+exp(logtheta)*H(beta))))
eval(%,[logtheta=1.5, H(beta)=2.0])

-- marginal survival (log-log link)
margS := log(-log((1+exp(logtheta)*H(beta))^(-exp(-logtheta))))
rule3 D(margS,beta)
rule3 D(margS,logtheta)