tests/startvals2.R

library(bbmle)

## fir data from emdbook package ...
firdata <- structure(list(TOTCONES = c(19, 42, 40, 68, 5, 0, 21, 114, 37, 
92, 84, 102, 98, 63, 9, 31, 35, 216, 27, 297, 36, 127, 23, 46, 
27, 66, 11, 20, 141, 3, 22, 39, 96, 206.5, 40, 231, 63.5, 202, 
54, 32, 107.5, 142.5, 82, 65, 153, 123, 131, 43, 98, 37, 34, 
10, 65, 35, 50, 19, 73, 33, 61, 9, 146, 0, 44, 42, 0, 61, 17, 
53, 27, 0, 74, 36, 28, 56, 46, 0, 15, 26, 46, 15, 105, 0, 62, 
24, 25, 41, 138, 77, 227.7, 28, 45, 57, 109, 28, 17, 91, 69, 
87, 10, 65, 50, 27, 30, 86, 119, 22, 8, 54, 104, 14, 16, 5, 53, 
40, 32, 114, 39, 37, 111, 226, 156, 42, 86, 94, 54, 1, 14, 44, 
108, 116.5, 14, 73, 3, 16, 87, 61, 48, 0, 17, 5, 88, 11, 133, 
121, 166, 171, 63, 23, 4, 51, 10, 14, 78, 47, 31, 42, 24, 42, 
55, 19, 63, 127, 9, 74, 120, 85, 51, 19, 131, 7, 23, 7, 9, 23, 
55, 48, 13, 2, 9, 3, 4, 16, 1, 88, 8, 27, 16, 184, 14, 22, 25, 
52, 2, 134, 81, 85, 3, 56, 17, 8, 10, 6, 69, 58, 1, 22, 3, 11, 
22, 2, 37, 8, 15, 61, 6, 18, 9, 109, 54, 4, 11, 30, 0, 0, 3, 
0, 16, 22, 9, 56, 17, 64, 38, 59, 37, 22, 41, 1, 22, 16, 17, 
4), DBH = c(9.4, 10.6, 7.7, 10.6, 8.7, 10.1, 8.1, 11.6, 10.1, 
13.3, 10, 13.4, 9.7, 7.4, 8.7, 8.6, 7.9, 14.2, 9.5, 15.9, 6, 
10.6, 7.3, 10.3, 8.4, 10.2, 13.8, 9.4, 8.1, 9.6, 7.3, 7.4, 10.3, 
13.4, 9.2, 13.9, 10.9, 17.4, 10.2, 8.2, 11.3, 16.1, 12.3, 8.3, 
12.4, 12.5, 11.3, 7.8, 11.6, 10, 7, 5.7, 7.7, 8.9, 8.5, 8.5, 
10.7, 10.2, 10.8, 9, 9.4, 7.6, 10.6, 10, 8, 7.4, 9.1, 6.7, 9.7, 
6.8, 8.6, 9.1, 6.3, 6.7, 10.9, 9.5, 9.9, 6.8, 9.8, 7.7, 12.1, 
8.2, 10, 9.6, 9.2, 8.2, 11.3, 11.6, 15.7, 9.1, 8.9, 8.7, 11, 
6.6, 7.1, 9, 12.4, 12.1, 7.5, 9, 8, 10.9, 9.2, 10.1, 12.1, 7, 
6.8, 8.6, 11.6, 6.6, 6.7, 6.8, 8.5, 7.8, 7.9, 9.8, 6.2, 6.7, 
15.4, 9.2, 12.9, 6.7, 9.6, 8.4, 8, 8.7, 6.7, 9.2, 9.5, 8, 5.5, 
8.5, 5.7, 5.6, 8, 6.5, 9.6, 6.1, 7.9, 5.9, 11, 8.2, 12.8, 12.8, 
12.5, 13.7, 11.8, 6.3, 6.3, 8.2, 6.2, 6.7, 9.8, 9.4, 6.7, 6, 
4.9, 9.6, 7.5, 8.4, 7.4, 9.9, 7.4, 9.5, 13.9, 6.9, 9.4, 7.4, 
12.8, 5.8, 7.2, 5.6, 6.9, 11.3, 9.6, 6.8, 6.9, 6.6, 4.8, 4.4, 
4.8, 8.5, 7, 8.7, 6.6, 8.6, 5.3, 10.4, 6.4, 5.4, 8.2, 5.5, 6.2, 
14.7, 10.5, 14.4, 5.8, 6.1, 6.2, 6.2, 7.2, 6, 10.6, 8.7, 7.5, 
7.3, 5.2, 6.9, 6.6, 6.7, 5.2, 6.9, 7.5, 9, 5.9, 6.5, 6.6, 9.8, 
4.7, 4.2, 4.8, 6.7, 6.5, 6.7, 5.9, 5.4, 6.9, 6.5, 6, 12, 7.5, 
6.4, 7.3, 7.3, 6.4, 7, 5.9, 9.1, 6.7, 4, 6.5, 4.7), WAVE_NON = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L), .Label = c("n", "w"), class = "factor"), logcones = c(2.99573227355399, 
3.76120011569356, 3.71357206670431, 4.23410650459726, 1.79175946922805, 
0, 3.09104245335832, 4.74493212836325, 3.63758615972639, 4.53259949315326, 
4.44265125649032, 4.63472898822964, 4.59511985013459, 4.15888308335967, 
2.30258509299405, 3.46573590279973, 3.58351893845611, 5.37989735354046, 
3.3322045101752, 5.6970934865054, 3.61091791264422, 4.85203026391962, 
3.17805383034795, 3.85014760171006, 3.3322045101752, 4.20469261939097, 
2.484906649788, 3.04452243772342, 4.95582705760126, 1.38629436111989, 
3.13549421592915, 3.68887945411394, 4.57471097850338, 5.33513133967075, 
3.71357206670431, 5.44673737166631, 4.16666522380173, 5.31320597904179, 
4.00733318523247, 3.49650756146648, 4.68675017298051, 4.96633503519968, 
4.4188406077966, 4.18965474202643, 5.03695260241363, 4.82028156560504, 
4.88280192258637, 3.78418963391826, 4.59511985013459, 3.63758615972639, 
3.55534806148941, 2.39789527279837, 4.18965474202643, 3.58351893845611, 
3.93182563272433, 2.99573227355399, 4.30406509320417, 3.52636052461616, 
4.12713438504509, 2.30258509299405, 4.99043258677874, 0, 3.80666248977032, 
3.76120011569356, 0, 4.12713438504509, 2.89037175789616, 3.98898404656427, 
3.3322045101752, 0, 4.31748811353631, 3.61091791264422, 3.36729582998647, 
4.04305126783455, 3.85014760171006, 0, 2.77258872223978, 3.29583686600433, 
3.85014760171006, 2.77258872223978, 4.66343909411207, 0, 4.14313472639153, 
3.2188758248682, 3.25809653802148, 3.73766961828337, 4.93447393313069, 
4.35670882668959, 5.43241110102874, 3.36729582998647, 3.8286413964891, 
4.06044301054642, 4.70048036579242, 3.36729582998647, 2.89037175789616, 
4.52178857704904, 4.24849524204936, 4.47733681447821, 2.39789527279837, 
4.18965474202643, 3.93182563272433, 3.3322045101752, 3.43398720448515, 
4.46590811865458, 4.78749174278205, 3.13549421592915, 2.19722457733622, 
4.00733318523247, 4.65396035015752, 2.70805020110221, 2.83321334405622, 
1.79175946922805, 3.98898404656427, 3.71357206670431, 3.49650756146648, 
4.74493212836325, 3.68887945411394, 3.63758615972639, 4.71849887129509, 
5.4249500174814, 5.05624580534831, 3.76120011569356, 4.46590811865458, 
4.55387689160054, 4.00733318523247, 0.693147180559945, 2.70805020110221, 
3.80666248977032, 4.69134788222914, 4.76643833358421, 2.70805020110221, 
4.30406509320417, 1.38629436111989, 2.83321334405622, 4.47733681447821, 
4.12713438504509, 3.89182029811063, 0, 2.89037175789616, 1.79175946922805, 
4.48863636973214, 2.484906649788, 4.89783979995091, 4.80402104473326, 
5.11799381241676, 5.14749447681345, 4.15888308335967, 3.17805383034795, 
1.6094379124341, 3.95124371858143, 2.39789527279837, 2.70805020110221, 
4.36944785246702, 3.87120101090789, 3.46573590279973, 3.76120011569356, 
3.2188758248682, 3.76120011569356, 4.02535169073515, 2.99573227355399, 
4.15888308335967, 4.85203026391962, 2.30258509299405, 4.31748811353631, 
4.79579054559674, 4.45434729625351, 3.95124371858143, 2.99573227355399, 
4.88280192258637, 2.07944154167984, 3.17805383034795, 2.07944154167984, 
2.30258509299405, 3.17805383034795, 4.02535169073515, 3.89182029811063, 
2.63905732961526, 1.09861228866811, 2.30258509299405, 1.38629436111989, 
1.6094379124341, 2.83321334405622, 0.693147180559945, 4.48863636973214, 
2.19722457733622, 3.3322045101752, 2.83321334405622, 5.22035582507832, 
2.70805020110221, 3.13549421592915, 3.25809653802148, 3.97029191355212, 
1.09861228866811, 4.90527477843843, 4.40671924726425, 4.45434729625351, 
1.38629436111989, 4.04305126783455, 2.89037175789616, 2.19722457733622, 
2.39789527279837, 1.94591014905531, 4.24849524204936, 4.07753744390572, 
0.693147180559945, 3.13549421592915, 1.38629436111989, 2.484906649788, 
3.13549421592915, 1.09861228866811, 3.63758615972639, 2.19722457733622, 
2.77258872223978, 4.12713438504509, 1.94591014905531, 2.94443897916644, 
2.30258509299405, 4.70048036579242, 4.00733318523247, 1.6094379124341, 
2.484906649788, 3.43398720448515, 0, 0, 1.38629436111989, 0, 
2.83321334405622, 3.13549421592915, 2.30258509299405, 4.04305126783455, 
2.89037175789616, 4.17438726989564, 3.66356164612965, 4.0943445622221, 
3.63758615972639, 3.13549421592915, 3.73766961828337, 0.693147180559945, 
3.13549421592915, 2.83321334405622, 2.89037175789616, 1.6094379124341
)), .Names = c("TOTCONES", "DBH", "WAVE_NON", "logcones"), row.names = c(1L, 
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 
16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 
29L, 30L, 31L, 32L, 33L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 
43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 
56L, 58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 
70L, 71L, 72L, 73L, 74L, 75L, 76L, 78L, 79L, 80L, 81L, 82L, 83L, 
84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L, 93L, 94L, 95L, 96L, 
97L, 98L, 99L, 100L, 101L, 102L, 103L, 104L, 105L, 106L, 107L, 
108L, 109L, 110L, 111L, 112L, 113L, 118L, 119L, 120L, 121L, 122L, 
123L, 124L, 126L, 127L, 128L, 129L, 130L, 131L, 132L, 133L, 134L, 
135L, 136L, 137L, 138L, 139L, 140L, 142L, 144L, 145L, 146L, 147L, 
148L, 149L, 150L, 151L, 154L, 155L, 157L, 159L, 160L, 168L, 169L, 
170L, 171L, 172L, 173L, 174L, 175L, 176L, 177L, 178L, 179L, 180L, 
181L, 184L, 185L, 186L, 187L, 189L, 190L, 193L, 198L, 247L, 272L, 
273L, 275L, 276L, 277L, 278L, 280L, 281L, 282L, 283L, 284L, 285L, 
286L, 287L, 288L, 289L, 290L, 291L, 292L, 293L, 294L, 295L, 296L, 
297L, 298L, 299L, 300L, 301L, 303L, 304L, 305L, 306L, 307L, 308L, 
309L, 310L, 311L, 313L, 314L, 315L, 316L, 319L, 320L, 321L, 322L, 
323L, 325L, 326L, 327L, 330L, 331L, 332L, 337L, 338L, 339L, 340L, 
341L, 342L, 343L, 344L, 345L, 346L, 347L, 348L, 349L, 350L, 351L, 
352L, 353L, 357L, 358L, 360L, 366L), na.action = structure(c(34L, 
57L, 77L, 114L, 115L, 116L, 117L, 125L, 141L, 143L, 152L, 153L, 
156L, 158L, 161L, 162L, 163L, 164L, 165L, 166L, 167L, 182L, 183L, 
188L, 191L, 192L, 194L, 195L, 196L, 197L, 199L, 200L, 201L, 202L, 
203L, 204L, 205L, 206L, 207L, 208L, 209L, 210L, 211L, 212L, 213L, 
214L, 215L, 216L, 217L, 218L, 219L, 220L, 221L, 222L, 223L, 224L, 
225L, 226L, 227L, 228L, 229L, 230L, 231L, 232L, 233L, 234L, 235L, 
236L, 237L, 238L, 239L, 240L, 241L, 242L, 243L, 244L, 245L, 246L, 
248L, 249L, 250L, 251L, 252L, 253L, 254L, 255L, 256L, 257L, 258L, 
259L, 260L, 261L, 262L, 263L, 264L, 265L, 266L, 267L, 268L, 269L, 
270L, 271L, 274L, 279L, 302L, 312L, 317L, 318L, 324L, 328L, 329L, 
333L, 334L, 335L, 336L, 354L, 355L, 356L, 359L, 361L, 362L, 363L, 
364L, 365L, 367L, 368L, 369L, 370L, 371L), .Names = c("34", "57", 
"77", "114", "115", "116", "117", "125", "141", "143", "152", 
"153", "156", "158", "161", "162", "163", "164", "165", "166", 
"167", "182", "183", "188", "191", "192", "194", "195", "196", 
"197", "199", "200", "201", "202", "203", "204", "205", "206", 
"207", "208", "209", "210", "211", "212", "213", "214", "215", 
"216", "217", "218", "219", "220", "221", "222", "223", "224", 
"225", "226", "227", "228", "229", "230", "231", "232", "233", 
"234", "235", "236", "237", "238", "239", "240", "241", "242", 
"243", "244", "245", "246", "248", "249", "250", "251", "252", 
"253", "254", "255", "256", "257", "258", "259", "260", "261", 
"262", "263", "264", "265", "266", "267", "268", "269", "270", 
"271", "274", "279", "302", "312", "317", "318", "324", "328", 
"329", "333", "334", "335", "336", "354", "355", "356", "359", 
"361", "362", "363", "364", "365", "367", "368", "369", "370", 
"371"), class = "omit"), class = "data.frame")


m1 <- mle2(logcones ~ dnorm(i + slope*log(DBH), sd),
           parameters= list(i ~ WAVE_NON-1, slope ~ WAVE_NON-1),
           data = firdata,
           start = list(i=c(-2,-2),slope=c(2.5,2.5),sd=1))

ancovafun = function(i1,i2,slope1,slope2,sigma) {
  int = c(i1,i2)[WAVE_NON]
  slope = c(slope1,slope2)[WAVE_NON]
  Y.pred =  int+ slope*log(DBH)
  r <- -sum(dnorm(logcones,mean=Y.pred,sd=sigma,log=TRUE))
  ## cat(i1,i2,slope1,slope2,sigma,r,"\n")
  r
}
m2 <- mle2(ancovafun,start=list(i1=-2,i2=-2,slope1=2.5,slope2=2.5,sigma=1),
  data=firdata)


m3 <- mle2(logcones ~ dnorm(mu, sd),
     parameters= list(mu ~ WAVE_NON*log(DBH)),
     data = firdata,
     start = list(mu=1,sd=1))

stopifnot(all.equal(AIC(m1),AIC(m2),AIC(m3)))

## m4 <- mle2(logcones ~ dnorm(i + slope*log(DBH), sd),
##           parameters= list(i ~ WAVE_NON-1, slope ~ WAVE_NON-1),
##           data = firdata,
##           start = c(-2,-2,2.5,2.5,sd=1))
bbolker/bbmle documentation built on Dec. 12, 2023, 9:07 a.m.