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