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]]