Dear R users,
I try to write an rcode for a paranetric weibull model with unobserved
heterogeneity.
The data I have are giuven below.
( The data comes from Keiding and Klein artricle published in
the Statistics in Medicine 1996) .
While the simple weibull model runs perfectly ( I get -2*loglikelihood
1093,905), the
weibull frailty model does not run well.
The code I write is :
likelihood.weibul <- function(p) {
cumhaz <- exp(x1*p[4])*(timeto^p[2])/p[1]
cumhaz <- sum(cumhaz)
lnhaz <- (x1*p[4]+log((p[2]*timeto^(p[2]-1))/p[1]))*stat
lik <- p[3]*log(p[3])-log(gamma(p[3]))+sum(lnhaz)+
log(gamma(di+p[3]))-(p[3]+di)*log(p[3]+cumhaz)
-2*lik
}
di <-sum(stat)
initial<-c(500000,2.8,.5,1.8)
# initial 1/exp(constant)
t<-nlm(likelihood.weibul,initial,print.level=1,hessian=T)
Best Regards,
Denis Lalountas
University of Patras
insem<-read.table("Melanoma.txt",header=T)
id <- insem$id
timeto<-insem$time
stat<-insem$status
x1<-insem$ulcer
x2 <-insem$thickness
x2 <- log(x2)
id time status sex age year thickness ulcer
1 1 10 0 1 76 1972 6.76 1
2 2 30 0 1 56 1968 0.65 0
3 3 35 0 1 41 1977 1.34 0
4 4 99 0 0 71 1968 2.90 0
5 5 185 1 1 52 1965 12.08 1
6 6 204 1 1 28 1971 4.84 1
7 7 210 1 1 77 1972 5.16 1
8 8 232 0 0 60 1974 3.22 1
9 9 232 1 1 49 1968 12.88 1
10 10 279 1 0 68 1971 7.41 1
11 11 295 1 0 53 1969 4.19 1
12 12 355 0 0 64 1972 0.16 1
13 13 386 1 0 68 1965 3.87 1
14 14 426 1 1 63 1970 4.84 1
15 15 469 1 0 14 1969 2.42 1
16 16 493 0 1 72 1971 12.56 1
17 17 529 1 1 46 1971 5.80 1
18 18 621 1 1 72 1972 7.06 1
19 19 629 1 1 95 1968 5.48 1
20 20 659 1 1 54 1972 7.73 1
21 21 667 1 0 89 1968 13.85 1
22 22 718 1 1 25 1967 2.34 1
23 23 752 1 1 37 1973 4.19 1
24 24 779 1 1 43 1967 4.04 1
25 25 793 1 1 68 1970 4.84 1
26 26 817 1 0 67 1966 0.32 0
27 27 826 0 0 86 1965 8.54 1
28 28 833 1 0 56 1971 2.58 1
29 29 858 1 0 16 1967 3.56 0
30 30 869 1 0 42 1965 3.54 0
31 31 872 1 0 65 1968 0.97 0
32 32 967 1 1 52 1970 4.83 1
33 33 977 1 1 58 1967 1.62 1
34 34 982 1 0 60 1970 6.44 1
35 35 1041 1 1 68 1967 14.66 0
36 36 1055 1 0 75 1967 2.58 1
37 37 1062 1 1 19 1966 3.87 1
38 38 1075 1 1 66 1971 3.54 1
39 39 1156 1 0 56 1970 1.34 1
40 40 1228 1 1 46 1973 2.24 1
41 41 1252 1 0 58 1971 3.87 1
42 42 1271 1 0 74 1971 3.54 1
43 43 1312 1 0 65 1970 17.42 1
44 44 1427 0 1 64 1972 1.29 0
45 45 1435 1 1 27 1969 3.22 0
46 46 1499 0 1 73 1973 1.29 0
47 47 1506 1 1 56 1970 4.51 1
48 48 1508 0 1 63 1973 8.38 1
49 49 1510 0 0 69 1973 1.94 0
50 50 1512 0 0 77 1973 0.16 0
51 51 1516 1 1 80 1968 2.58 1
52 52 1525 0 0 76 1970 1.29 1
53 53 1542 0 0 65 1973 0.16 0
54 54 1548 1 0 61 1972 1.62 0
55 55 1557 0 0 26 1973 1.29 0
56 56 1560 1 0 57 1973 2.10 0
57 57 1563 0 0 45 1973 0.32 0
58 58 1584 1 1 31 1970 0.81 0
59 59 1605 0 0 36 1973 1.13 0
60 60 1621 1 0 46 1972 5.16 1
61 61 1627 0 0 43 1973 1.62 0
62 62 1634 0 0 68 1973 1.37 0
63 63 1641 0 1 57 1973 0.24 0
64 64 1641 0 0 57 1973 0.81 0
65 65 1648 0 0 55 1973 1.29 0
66 66 1652 0 0 58 1973 1.29 0
67 67 1654 0 1 20 1973 0.97 0
68 68 1654 0 0 67 1973 1.13 0
69 69 1667 1 0 44 1971 5.80 1
70 70 1678 0 0 59 1973 1.29 0
71 71 1685 0 0 32 1973 0.48 0
72 72 1690 1 1 83 1971 1.62 0
73 73 1710 0 0 55 1973 2.26 0
74 74 1710 0 1 15 1973 0.58 0
75 75 1726 1 0 58 1970 0.97 1
76 76 1745 0 0 47 1973 2.58 1
77 77 1762 0 0 54 1973 0.81 0
78 78 1779 0 1 55 1973 3.54 1
79 79 1787 0 1 38 1973 0.97 0
80 80 1787 0 0 41 1973 1.78 1
81 81 1793 0 0 56 1973 1.94 0
82 82 1804 0 0 48 1973 1.29 0
83 83 1812 0 1 44 1973 3.22 1
84 84 1836 0 0 70 1972 1.53 0
85 85 1839 0 0 40 1972 1.29 0
86 86 1839 0 1 53 1972 1.62 1
87 87 1854 0 0 65 1972 1.62 1
88 88 1856 0 1 54 1972 0.32 0
89 89 1860 0 1 71 1969 4.84 1
90 90 1864 0 0 49 1972 1.29 0
91 91 1899 0 0 55 1972 0.97 0
92 92 1914 0 0 69 1972 3.06 0
93 93 1919 0 1 83 1972 3.54 0
94 94 1920 0 1 60 1972 1.62 1
95 95 1927 0 1 40 1972 2.58 1
96 96 1933 1 0 77 1972 1.94 0
97 97 1942 0 0 35 1972 0.81 0
98 98 1955 0 0 46 1972 7.73 1
99 99 1956 0 0 34 1972 0.97 0
100 100 1958 0 0 69 1972 12.88 0
101 101 1963 0 0 60 1972 2.58 0
102 102 1970 0 1 84 1972 4.09 1
103 103 2005 0 0 66 1972 0.64 0
104 104 2007 0 1 56 1972 0.97 0
105 105 2011 0 0 75 1972 3.22 1
106 106 2024 0 0 36 1972 1.62 0
107 107 2028 0 1 52 1972 3.87 1
108 108 2038 0 0 58 1972 0.32 1
109 109 2056 0 0 39 1972 0.32 0
110 110 2059 0 1 68 1972 3.22 1
111 111 2061 1 1 71 1968 2.26 0
112 112 2062 1 0 52 1965 3.06 0
113 113 2075 0 1 55 1972 2.58 1
114 114 2085 0 0 66 1970 0.65 0
115 115 2102 0 1 35 1972 1.13 0
116 116 2103 1 1 44 1966 0.81 0
117 117 2104 0 0 72 1972 0.97 0
118 118 2108 1 0 58 1969 1.76 1
119 119 2112 0 0 54 1972 1.94 1
120 120 2150 0 0 33 1972 0.65 0
121 121 2156 0 0 45 1972 0.97 0
122 122 2165 0 1 62 1972 5.64 0
123 123 2209 0 0 72 1971 9.66 0
124 124 2227 0 0 51 1971 0.10 0
125 125 2227 0 1 77 1971 5.48 1
126 126 2256 1 0 43 1971 2.26 1
127 127 2264 0 0 65 1971 4.83 1
128 128 2339 0 0 63 1971 0.97 0
129 129 2361 0 1 60 1971 0.97 0
130 130 2387 0 0 50 1971 5.16 1
131 131 2388 1 1 40 1966 0.81 0
132 132 2403 0 0 67 1971 2.90 1
133 133 2426 0 0 69 1971 3.87 0
134 134 2426 0 0 74 1971 1.94 1
135 135 2431 0 0 49 1971 0.16 0
136 136 2460 0 0 47 1971 0.64 0
137 137 2467 1 0 42 1965 2.26 1
138 138 2492 0 0 54 1971 1.45 0
139 139 2493 0 1 72 1971 4.82 1
140 140 2521 0 0 45 1971 1.29 1
141 141 2542 0 1 67 1971 7.89 1
142 142 2559 0 0 48 1970 0.81 1
143 143 2565 1 1 34 1970 3.54 1
144 144 2570 0 0 44 1970 1.29 0
145 145 2660 0 0 31 1970 0.64 0
146 146 2666 0 0 42 1970 3.22 1
147 147 2676 0 0 24 1970 1.45 1
148 148 2738 0 0 58 1970 0.48 0
149 149 2782 1 1 78 1969 1.94 0
150 150 2787 0 1 62 1970 0.16 0
151 151 2984 0 1 70 1969 0.16 0
152 152 3032 0 0 35 1969 1.29 0
153 153 3040 0 0 61 1969 1.94 0
154 154 3042 1 0 54 1967 3.54 1
155 155 3067 0 0 29 1969 0.81 0
156 156 3079 0 1 64 1969 0.65 0
157 157 3101 0 1 47 1969 7.09 0
158 158 3144 0 1 62 1969 0.16 0
159 159 3152 0 0 32 1969 1.62 0
160 160 3154 0 1 49 1969 1.62 0
161 161 3180 0 0 25 1969 1.29 0
162 162 3182 0 1 49 1966 6.12 0
163 163 3185 0 0 64 1969 0.48 0
164 164 3199 0 0 36 1969 0.64 0
165 165 3228 0 0 58 1969 3.22 1
166 166 3229 0 0 37 1969 1.94 0
167 167 3278 0 1 54 1969 2.58 0
168 168 3297 0 0 61 1968 2.58 1
169 169 3328 0 1 31 1968 0.81 0
170 170 3330 0 1 61 1968 0.81 1
171 171 3338 1 0 60 1967 3.22 1
172 172 3383 0 0 43 1968 0.32 0
173 173 3384 0 0 68 1968 3.22 1
174 174 3385 0 0 4 1968 2.74 0
175 175 3388 0 1 60 1968 4.84 1
176 176 3402 0 1 50 1968 1.62 0
177 177 3441 0 0 20 1968 0.65 0
178 178 3458 0 0 54 1967 1.45 0
179 179 3459 0 0 29 1968 0.65 0
180 180 3459 0 1 56 1968 1.29 1
181 181 3476 0 0 60 1968 1.62 0
182 182 3523 0 0 46 1968 3.54 0
183 183 3667 0 0 42 1967 3.22 0
184 184 3695 0 0 34 1967 0.65 0
185 185 3695 0 0 56 1967 1.03 0
186 186 3776 0 1 12 1967 7.09 1
187 187 3776 0 0 21 1967 1.29 1
188 188 3830 0 1 46 1967 0.65 0
189 189 3856 0 0 49 1967 1.78 0
190 190 3872 0 0 35 1967 12.24 1
191 191 3909 0 1 42 1967 8.06 1
192 192 3968 0 0 47 1967 0.81 0
193 193 4001 0 0 69 1967 2.10 0
194 194 4103 0 0 52 1966 3.87 0
195 195 4119 0 1 52 1966 0.65 0
196 196 4124 0 0 30 1966 1.94 1
197 197 4207 0 1 22 1966 0.65 0
198 198 4310 0 1 55 1966 2.10 0
199 199 4390 0 0 26 1965 1.94 1
200 200 4479 0 0 19 1965 1.13 1
201 201 4492 0 1 29 1965 7.06 1
202 202 4668 0 0 40 1965 6.12 0
203 203 4688 0 0 42 1965 0.48 0
204 204 4926 0 0 50 1964 2.26 0
205 205 5565 0 0 41 1962 2.90 0
---------------------------------
Boardwalk for $500? In 2007? Ha!
[[alternative HTML version deleted]]
