[Lme4-commits] r1613 - pkg/lme4Eigen/tests

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Feb 23 19:39:03 CET 2012


Author: bbolker
Date: 2012-02-23 19:39:03 +0100 (Thu, 23 Feb 2012)
New Revision: 1613

Added:
   pkg/lme4Eigen/tests/Elston2001_tickdata.txt
   pkg/lme4Eigen/tests/elston.R
Log:

  add elston data and example
  (should eventually move data to data/ directory, document, etc. --
or move to mlmRev)



Added: pkg/lme4Eigen/tests/Elston2001_tickdata.txt
===================================================================
--- pkg/lme4Eigen/tests/Elston2001_tickdata.txt	                        (rev 0)
+++ pkg/lme4Eigen/tests/Elston2001_tickdata.txt	2012-02-23 18:39:03 UTC (rev 1613)
@@ -0,0 +1,405 @@
+	INDEX	 TICKS	 BROOD	 HEIGHT	 YEAR    LOCATION
+	1	 0	 501	 465	 95	 32
+	2	 0	 501	 465	 95	 32
+	3	 0	 502	 472	 95	 36
+	4	 0	 503	 475	 95	 37
+	5	 0	 503	 475	 95	 37
+	6	 3	 503	 475	 95	 37
+	7	 2	 503	 475	 95	 37
+	8	 0	 504	 488	 95	 44
+	9	 0	 504	 488	 95	 44
+	10	 2	 504	 488	 95	 44
+	11	 0	 505	 492	 95	 47
+	12	 0	 505	 492	 95	 47
+	13	 0	 505	 492	 95	 47
+	14	 0	 506	 490	 95	 45
+	15	 0	 506	 490	 95	 45
+	16	 0	 506	 490	 95	 45
+	17	 0	 507	 464	 95	 31
+	18	 0	 507	 464	 95	 31
+	19	 0	 507	 464	 95	 31
+	20	 1	 507	 464	 95	 31
+	21	 2	 507	 464	 95	 31
+	22	 1	 509	 457	 95	 28
+	23	 0	 510	 457	 95	 28
+	24	 0	 511	 457	 95	 28
+	25	 5	 511	 457	 95	 28
+	26	 8	 512	 451	 95	 26
+	27	 3	 512	 451	 95	 26
+	28	 4	 512	 451	 95	 26
+	29	 7	 513	 437	 95	 17
+	30	 0	 513	 437	 95	 17
+	31	 4	 513	 437	 95	 17
+	32	 4	 514	 430	 95	 14
+	33	 1	 514	 430	 95	 14
+	34	 0	 514	 430	 95	 14
+	35	 0	 514	 430	 95	 14
+	36	 3	 514	 430	 95	 14
+	37	 1	 514	 430	 95	 14
+	38	 6	 515	 427	 95	 13
+	39	 0	 515	 427	 95	 13
+	40	 1	 515	 427	 95	 13
+	41	 0	 515	 427	 95	 13
+	42	 2	 516	 419	 95	 7
+	43	 7	 516	 419	 95	 7
+	44	 31	 516	 419	 95	 7
+	45	 34	 517	 411	 95	 4
+	46	 17	 517	 411	 95	 4
+	47	 16	 517	 411	 95	 4
+	48	 66	 518	 419	 95	 7
+	49	 49	 518	 419	 95	 7
+	50	 82	 518	 419	 95	 7
+	51	 85	 518	 419	 95	 7
+	52	 64	 518	 419	 95	 7
+	53	 11	 519	 424	 95	 11
+	54	 14	 519	 424	 95	 11
+	55	 4	 519	 424	 95	 11
+	56	 10	 519	 424	 95	 11
+	57	 3	 520	 427	 95	 13
+	58	 15	 520	 427	 95	 13
+	59	 8	 520	 427	 95	 13
+	60	 9	 521	 422	 95	 9
+	61	 11	 521	 422	 95	 9
+	62	 7	 521	 422	 95	 9
+	63	 13	 521	 422	 95	 9
+	64	 3	 522	 503	 95	 53
+	65	 0	 523	 496	 95	 51
+	66	 0	 523	 496	 95	 51
+	67	 1	 523	 496	 95	 51
+	68	 1	 523	 496	 95	 51
+	69	 0	 525	 507	 95	 54
+	70	 0	 525	 507	 95	 54
+	71	 0	 525	 507	 95	 54
+	72	 0	 525	 507	 95	 54
+	73	 0	 526	 496	 95	 51
+	74	 0	 526	 496	 95	 51
+	75	 0	 526	 496	 95	 51
+	76	 1	 526	 496	 95	 51
+	77	 1	 526	 496	 95	 51
+	78	 0	 526	 496	 95	 51
+	79	 2	 528	 466	 95	 33
+	80	 0	 528	 466	 95	 33
+	81	 3	 528	 466	 95	 33
+	82	 1	 531	 488	 95	 44
+	83	 7	 533	 442	 95	 19
+	84	 2	 533	 442	 95	 19
+	85	 16	 533	 442	 95	 19
+	86	 12	 533	 442	 95	 19
+	87	 0	 533	 442	 95	 19
+	88	 1	 535	 442	 95	 19
+	89	 0	 537	 533	 95	 63
+	90	 0	 537	 533	 95	 63
+	91	 1	 537	 533	 95	 63
+	92	 0	 537	 533	 95	 63
+	93	 0	 537	 533	 95	 63
+	94	 1	 537	 533	 95	 63
+	95	 0	 537	 533	 95	 63
+	96	 0	 538	 533	 95	 63
+	97	 0	 539	 515	 95	 59
+	98	 0	 539	 515	 95	 59
+	99	 0	 539	 515	 95	 59
+	100	 0	 539	 515	 95	 59
+	101	 5	 540	 518	 95	 60
+	102	 2	 540	 518	 95	 60
+	103	 2	 542	 493	 95	 48
+	104	 1	 542	 493	 95	 48
+	105	 1	 542	 493	 95	 48
+	106	 0	 548	 468	 95	 34
+	107	 0	 548	 468	 95	 34
+	108	 1	 548	 468	 95	 34
+	109	 1	 549	 476	 95	 38
+	110	 1	 549	 476	 95	 38
+	111	 0	 549	 476	 95	 38
+	112	 0	 549	 476	 95	 38
+	113	 5	 550	 446	 95	 22
+	114	 3	 550	 446	 95	 22
+	115	 2	 553	 460	 95	 30
+	116	 2	 559	 525	 95	 62
+	117	 1	 559	 525	 95	 62
+	118	 1	 601	 410	 96	 3
+	119	 0	 601	 410	 96	 3
+	120	 2	 601	 410	 96	 3
+	121	 5	 601	 410	 96	 3
+	122	 1	 601	 410	 96	 3
+	123	 2	 601	 410	 96	 3
+	124	 2	 601	 410	 96	 3
+	125	 3	 602	 417	 96	 6
+	126	 14	 602	 417	 96	 6
+	127	 11	 602	 417	 96	 6
+	128	 9	 602	 417	 96	 6
+	129	 4	 602	 417	 96	 6
+	130	 10	 602	 417	 96	 6
+	131	 33	 602	 417	 96	 6
+	132	 19	 602	 417	 96	 6
+	133	 16	 602	 417	 96	 6
+	134	 16	 603	 430	 96	 14
+	135	 13	 603	 430	 96	 14
+	136	 11	 603	 430	 96	 14
+	137	 7	 603	 430	 96	 14
+	138	 4	 603	 430	 96	 14
+	139	 11	 603	 430	 96	 14
+	140	 1	 604	 456	 96	 27
+	141	 1	 604	 456	 96	 27
+	142	 4	 604	 456	 96	 27
+	143	 6	 605	 457	 96	 28
+	144	 2	 605	 457	 96	 28
+	145	 7	 605	 457	 96	 28
+	146	 8	 605	 457	 96	 28
+	147	 14	 605	 457	 96	 28
+	148	 6	 606	 430	 96	 14
+	149	 13	 606	 430	 96	 14
+	150	 5	 606	 430	 96	 14
+	151	 8	 606	 430	 96	 14
+	152	 13	 606	 430	 96	 14
+	153	 17	 606	 430	 96	 14
+	154	 5	 606	 430	 96	 14
+	155	 1	 606	 430	 96	 14
+	156	 1	 606	 430	 96	 14
+	157	 2	 606	 430	 96	 14
+	158	 7	 607	 423	 96	 10
+	159	 7	 608	 421	 96	 8
+	160	 11	 608	 421	 96	 8
+	161	 1	 609	 525	 96	 62
+	162	 0	 609	 525	 96	 62
+	163	 5	 610	 509	 96	 55
+	164	 4	 610	 509	 96	 55
+	165	 0	 611	 499	 96	 52
+	166	 0	 611	 499	 96	 52
+	167	 0	 611	 499	 96	 52
+	168	 7	 612	 503	 96	 53
+	169	 5	 612	 503	 96	 53
+	170	 3	 612	 503	 96	 53
+	171	 1	 612	 503	 96	 53
+	172	 6	 612	 503	 96	 53
+	173	 1	 614	 492	 96	 47
+	174	 2	 614	 492	 96	 47
+	175	 14	 615	 491	 96	 46
+	176	 5	 615	 491	 96	 46
+	177	 27	 615	 491	 96	 46
+	178	 1	 616	 475	 96	 37
+	179	 2	 616	 475	 96	 37
+	180	 3	 616	 475	 96	 37
+	181	 0	 616	 475	 96	 37
+	182	 1	 617	 479	 96	 40
+	183	 0	 617	 479	 96	 40
+	184	 5	 617	 479	 96	 40
+	185	 5	 617	 479	 96	 40
+	186	 8	 617	 479	 96	 40
+	187	 21	 617	 479	 96	 40
+	188	 15	 618	 472	 96	 36
+	189	 15	 618	 472	 96	 36
+	190	 6	 618	 472	 96	 36
+	191	 19	 618	 472	 96	 36
+	192	 14	 618	 472	 96	 36
+	193	 1	 621	 485	 96	 42
+	194	 1	 621	 485	 96	 42
+	195	 3	 621	 485	 96	 42
+	196	 2	 621	 485	 96	 42
+	197	 3	 621	 485	 96	 42
+	198	 2	 623	 495	 96	 50
+	199	 5	 623	 495	 96	 50
+	200	 0	 624	 472	 96	 36
+	201	 6	 624	 472	 96	 36
+	202	 3	 624	 472	 96	 36
+	203	 1	 625	 458	 96	 29
+	204	 0	 625	 458	 96	 29
+	205	 1	 625	 458	 96	 29
+	206	 6	 625	 458	 96	 29
+	207	 1	 625	 458	 96	 29
+	208	 85	 626	 449	 96	 24
+	209	 45	 626	 449	 96	 24
+	210	 68	 626	 449	 96	 24
+	211	 84	 626	 449	 96	 24
+	212	 50	 626	 449	 96	 24
+	213	 13	 628	 442	 96	 19
+	214	 1	 628	 442	 96	 19
+	215	 19	 629	 448	 96	 23
+	216	 26	 629	 448	 96	 23
+	217	 9	 629	 448	 96	 23
+	218	 2	 629	 448	 96	 23
+	219	 4	 629	 448	 96	 23
+	220	 3	 629	 448	 96	 23
+	221	 22	 630	 448	 96	 23
+	222	 32	 630	 448	 96	 23
+	223	 5	 631	 403	 96	 1
+	224	 21	 631	 403	 96	 1
+	225	 26	 631	 403	 96	 1
+	226	 13	 631	 403	 96	 1
+	227	 23	 631	 403	 96	 1
+	228	 42	 632	 411	 96	 4
+	229	 38	 632	 411	 96	 4
+	230	 61	 632	 411	 96	 4
+	231	 79	 632	 411	 96	 4
+	232	 39	 632	 411	 96	 4
+	233	 41	 632	 411	 96	 4
+	234	 15	 634	 415	 96	 5
+	235	 23	 634	 415	 96	 5
+	236	 14	 634	 415	 96	 5
+	237	 7	 635	 427	 96	 13
+	238	 24	 636	 424	 96	 11
+	239	 3	 638	 525	 96	 62
+	240	 1	 638	 525	 96	 62
+	241	 2	 640	 521	 96	 61
+	242	 1	 640	 521	 96	 61
+	243	 0	 640	 521	 96	 61
+	244	 3	 641	 518	 96	 60
+	245	 8	 641	 518	 96	 60
+	246	 1	 642	 495	 96	 50
+	247	 2	 642	 495	 96	 50
+	248	 0	 642	 495	 96	 50
+	249	 8	 643	 495	 96	 50
+	250	 3	 643	 495	 96	 50
+	251	 14	 643	 495	 96	 50
+	252	 16	 643	 495	 96	 50
+	253	 18	 643	 495	 96	 50
+	254	 11	 643	 495	 96	 50
+	255	 13	 643	 495	 96	 50
+	256	 6	 645	 460	 96	 30
+	257	 7	 645	 460	 96	 30
+	258	 10	 645	 460	 96	 30
+	259	 5	 647	 442	 96	 19
+	260	 7	 647	 442	 96	 19
+	261	 25	 648	 443	 96	 20
+	262	 11	 648	 443	 96	 20
+	263	 6	 648	 443	 96	 20
+	264	 4	 648	 443	 96	 20
+	265	 7	 648	 443	 96	 20
+	266	 4	 650	 425	 96	 12
+	267	 6	 650	 425	 96	 12
+	268	 2	 650	 425	 96	 12
+	269	 5	 651	 439	 96	 18
+	270	 3	 651	 439	 96	 18
+	271	 7	 651	 439	 96	 18
+	272	 3	 652	 444	 96	 21
+	273	 1	 701	 450	 97	 25
+	274	 4	 701	 450	 97	 25
+	275	 4	 701	 450	 97	 25
+	276	 2	 701	 450	 97	 25
+	277	 5	 701	 450	 97	 25
+	278	 3	 702	 446	 97	 22
+	279	 0	 702	 446	 97	 22
+	280	 3	 702	 446	 97	 22
+	281	 1	 702	 446	 97	 22
+	282	 2	 702	 446	 97	 22
+	283	 3	 702	 446	 97	 22
+	284	 1	 704	 472	 97	 36
+	285	 0	 704	 472	 97	 36
+	286	 4	 704	 472	 97	 36
+	287	 0	 704	 472	 97	 36
+	288	 0	 704	 472	 97	 36
+	289	 0	 705	 472	 97	 36
+	290	 0	 706	 460	 97	 30
+	291	 0	 706	 460	 97	 30
+	292	 0	 706	 460	 97	 30
+	293	 1	 708	 442	 97	 19
+	294	 3	 708	 442	 97	 19
+	295	 4	 708	 442	 97	 19
+	296	 0	 708	 442	 97	 19
+	297	 4	 708	 442	 97	 19
+	298	 2	 708	 442	 97	 19
+	299	 0	 709	 525	 97	 62
+	300	 0	 709	 525	 97	 62
+	301	 1	 709	 525	 97	 62
+	302	 0	 710	 533	 97	 63
+	303	 1	 710	 533	 97	 63
+	304	 2	 710	 533	 97	 63
+	305	 0	 710	 533	 97	 63
+	306	 0	 710	 533	 97	 63
+	307	 0	 710	 533	 97	 63
+	308	 1	 711	 513	 97	 57
+	309	 0	 711	 513	 97	 57
+	310	 0	 711	 513	 97	 57
+	311	 1	 711	 513	 97	 57
+	312	 1	 711	 513	 97	 57
+	313	 0	 711	 513	 97	 57
+	314	 0	 711	 513	 97	 57
+	315	 0	 712	 514	 97	 58
+	316	 0	 712	 514	 97	 58
+	317	 0	 713	 511	 97	 56
+	318	 0	 713	 511	 97	 56
+	319	 1	 713	 511	 97	 56
+	320	 0	 713	 511	 97	 56
+	321	 0	 713	 511	 97	 56
+	322	 0	 714	 511	 97	 56
+	323	 0	 714	 511	 97	 56
+	324	 1	 714	 511	 97	 56
+	325	 0	 714	 511	 97	 56
+	326	 0	 714	 511	 97	 56
+	327	 0	 714	 511	 97	 56
+	328	 0	 715	 496	 97	 51
+	329	 0	 715	 496	 97	 51
+	330	 0	 715	 496	 97	 51
+	331	 0	 715	 496	 97	 51
+	332	 1	 716	 494	 97	 49
+	333	 0	 716	 494	 97	 49
+	334	 0	 717	 494	 97	 49
+	335	 0	 717	 494	 97	 49
+	336	 0	 717	 494	 97	 49
+	337	 2	 718	 411	 97	 4
+	338	 4	 718	 411	 97	 4
+	339	 4	 718	 411	 97	 4
+	340	 1	 718	 411	 97	 4
+	341	 2	 718	 411	 97	 4
+	342	 2	 719	 411	 97	 4
+	343	 1	 719	 411	 97	 4
+	344	 4	 719	 411	 97	 4
+	345	 1	 719	 411	 97	 4
+	346	 3	 719	 411	 97	 4
+	347	 3	 719	 411	 97	 4
+	348	 2	 720	 423	 97	 10
+	349	 0	 720	 423	 97	 10
+	350	 3	 721	 424	 97	 11
+	351	 2	 722	 403	 97	 1
+	352	 0	 722	 403	 97	 1
+	353	 1	 722	 403	 97	 1
+	354	 0	 723	 409	 97	 2
+	355	 3	 723	 409	 97	 2
+	356	 1	 723	 409	 97	 2
+	357	 10	 724	 434	 97	 16
+	358	 4	 724	 434	 97	 16
+	359	 1	 724	 434	 97	 16
+	360	 2	 725	 477	 97	 39
+	361	 0	 725	 477	 97	 39
+	362	 1	 725	 477	 97	 39
+	363	 0	 727	 472	 97	 36
+	364	 0	 728	 468	 97	 34
+	365	 0	 728	 468	 97	 34
+	366	 0	 729	 470	 97	 35
+	367	 0	 730	 486	 97	 43
+	368	 0	 731	 495	 97	 50
+	369	 0	 731	 495	 97	 50
+	370	 0	 731	 495	 97	 50
+	371	 0	 731	 495	 97	 50
+	372	 0	 731	 495	 97	 50
+	373	 0	 732	 483	 97	 41
+	374	 0	 732	 483	 97	 41
+	375	 2	 732	 483	 97	 41
+	376	 2	 733	 442	 97	 19
+	377	 2	 733	 442	 97	 19
+	378	 0	 734	 457	 97	 28
+	379	 4	 736	 457	 97	 28
+	380	 5	 736	 457	 97	 28
+	381	 2	 737	 457	 97	 28
+	382	 3	 737	 457	 97	 28
+	383	 2	 737	 457	 97	 28
+	384	 2	 737	 457	 97	 28
+	385	 1	 737	 457	 97	 28
+	386	 2	 737	 457	 97	 28
+	387	 0	 737	 457	 97	 28
+	388	 2	 738	 464	 97	 31
+	389	 0	 738	 464	 97	 31
+	390	 1	 738	 464	 97	 31
+	391	 0	 739	 433	 97	 15
+	392	 3	 739	 433	 97	 15
+	393	 1	 739	 433	 97	 15
+	394	 0	 739	 433	 97	 15
+	395	 0	 740	 442	 97	 19
+	396	 0	 740	 442	 97	 19
+	397	 0	 741	 433	 97	 15
+	398	 1	 741	 433	 97	 15
+	399	 0	 741	 433	 97	 15
+	400	 0	 742	 430	 97	 14
+	401	 0	 742	 430	 97	 14
+	402	 2	 743	 450	 97	 25
+	403	 0	 743	 450	 97	 25
+ 

Added: pkg/lme4Eigen/tests/elston.R
===================================================================
--- pkg/lme4Eigen/tests/elston.R	                        (rev 0)
+++ pkg/lme4Eigen/tests/elston.R	2012-02-23 18:39:03 UTC (rev 1613)
@@ -0,0 +1,61 @@
+tickdata <- read.table("Elston2001_tickdata.txt",header=TRUE,
+  colClasses=c("factor","numeric","factor","numeric","factor","factor"))
+
+tickdata <- transform(tickdata,cHEIGHT=scale(HEIGHT,scale=FALSE))
+
+form <- TICKS~YEAR+HEIGHT+(1|BROOD)+(1|INDEX)+(1|LOCATION)
+## fit with lme4
+library(lme4)
+t1 <- system.time(full_mod1  <- glmer(form, family="poisson",data=tickdata))
+c1 <- c(fixef(full_mod1),unlist(VarCorr(full_mod1)), logLik=logLik(full_mod1),time=t1["elapsed"])
+allcoefs1 <- c(unlist(full_mod1 at ST),fixef(full_mod1))
+detach("package:lme4")
+
+## fit with lme4Eigen
+library(lme4Eigen)
+t2 <- system.time(full_mod2  <- glmer(form, family="poisson",data=tickdata))
+c2 <- c(fixef(full_mod2),unlist(VarCorr(full_mod2)), logLik=logLik(full_mod2),time=t2["elapsed"])
+
+allcoefs <- function(x) c(getME(x,"theta"),getME(x,"beta"))
+
+## deviance function
+mm <- glmer(form, family="poisson",data=tickdata,devFunOnly=TRUE,tolPwrss=1e-12,verbose=4,
+            compDev=FALSE)
+mm(allcoefs1)
+## works with compDev=FALSE, fails with compDev=TRUE
+
+## refit
+full_mod3 <- refit(full_mod2,tickdata$TICKS)
+
+fn <- "elston_fits.RData"
+## FIXME::: some possibility of differing results? 1780.
+## what changed ???
+if (!file.exists(fn)) {
+  tvec <- seq(-13,-7,by=0.1)
+  dmat <- matrix(nrow=length(tvec),ncol=9,  
+                 dimnames=list(NULL,c("deviance","time_elapsed",
+                   paste("theta",1:3,sep=""),paste("beta",1:4,sep=""))))
+  for (i in seq_along(tvec)) {
+    tt <- system.time(gg <- glmer(form,family="poisson",data=tickdata,tolPwrss=10^tvec[i]))["elapsed"]
+    cat(i,tvec[i],deviance(gg),"\n")
+    dmat[i,] <- c(deviance(gg),tt,allcoefs(gg))
+  }
+  dmat <- data.frame(logtol=tvec,dmat)
+  save("dmat",file=fn)
+} else load(fn)
+
+newdev <- apply(dmat[,-(1:3)],1,mm)
+library(ggplot2)
+library(reshape)
+qplot(logtol,value,data=melt(dmat,id.var="logtol"),geom=c("line"))+
+  facet_wrap(~variable,scale="free")+
+  geom_line(data=data.frame(logtol=dmat$logtol,value=newdev,variable="deviance"),
+            colour="blue")+
+  geom_hline(data=data.frame(variable="deviance",val=mm(allcoefs1)),
+             aes(yintercept=val),colour="red")+theme_bw()+
+  geom_vline(data=data.frame(variable="deviance",val=-10),
+             aes(xintercept=val),colour="gray")
+
+detach("package:lme4Eigen")
+cbind(c1,c2)
+



More information about the Lme4-commits mailing list