В прошлом я использовал asreml-r для учета пространственной автокорреляции в полевых испытаниях в сельском хозяйстве, которые были построены по схеме «ряд и диапазон».Относительно легко использовать пакет asreml для указания пространственной модели (то есть rcov = ~ at (LOCATION): ar1 (ROW): ar1 (RANGE))
К сожалению, asreml-r дорог и труден дляучить.Моя исследовательская группа также предпочитает полагаться на nlme и lmer для большинства своих аналитических потребностей.Поэтому они не хотят платить за asreml-r или рассмотреть возможность использования.
Несколько лет назад был опубликован вопрос с вопросом, доступна ли альтернатива asreml-r с открытым исходным кодом, которую можно использовать для построения двумерной пространственной модели со структурой ошибок в обоих направлениях.,В то время казалось, что все согласны с тем, что делать это либо в lmer, либо в nlme было не просто.
Проведя несколько часов в поисках, я не совсем понимаю, был ли достигнут какой-либо прогресс в этой области.обращаясь к этому.Кто-нибудь может отослать меня к недавней дискуссии относительно этого типа анализа?Или они могут дать совет о том, как построить модели со смешанными эффектами, которые учитывают пространственную корреляцию в nlme или lmer?
Обратите внимание, что ни я, ни другие члены нашей группы не являются в точности статистиками или высококвалифицированными программистами.Также непрактично заключать контракт с внешней группой для анализа наших данных.Мы просто хотим применять лучшие методы, которые мы можем использовать для регулярного ежегодного анализа данных.
Пример анализируемых данных:
my.data <- structure(list(ENTRY = structure(c(23L, 23L, 23L, 40L, 12L, 8L,
1L, 15L, 30L, 1L, 24L, 8L, 1L, 8L, 30L, 33L, 12L, 38L, 41L, 36L,
43L, 32L, 44L, 31L, 26L, 11L, 13L, 34L, 5L, 22L, 4L, 14L, 11L,
20L, 25L, 11L, 21L, 43L, 44L, 4L, 42L, 45L, 42L, 41L, 42L, 4L,
44L, 20L, 40L, 29L, 29L, 24L, 2L, 3L, 28L, 24L, 34L, 27L, 41L,
28L, 29L, 5L, 3L, 25L, 14L, 20L, 15L, 21L, 31L, 22L, 40L, 21L,
6L, 38L, 43L, 12L, 6L, 14L, 5L, 3L, 30L, 45L, 31L, 7L, 9L, 39L,
22L, 15L, 26L, 28L, 34L, 10L, 25L, 27L, 16L, 45L, 10L, 18L, 32L,
10L, 6L, 18L, 33L, 16L, 37L, 9L, 32L, 38L, 39L, 2L, 2L, 39L,
36L, 36L, 7L, 27L, 7L, 26L, 17L, 9L, 33L, 13L, 17L, 17L, 35L,
37L, 37L, 18L, 16L, 19L, 13L, 19L, 35L, 19L, 35L, 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, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L,
44L, 45L, 52L, 54L, 52L, 54L, 49L, 51L, 50L, 54L, 49L, 46L, 51L,
50L, 53L, 49L, 50L, 51L, 53L, 52L, 53L, 48L, 47L, 46L, 46L, 47L,
48L, 48L, 47L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L), .Label = c("20",
"112", "1478", "1495", "1521", "1522", "1590", "1608", "1657",
"1658", "1660", "1667", "1680", "1688", "1723", "1728", "1730",
"1731", "1743", "1745", "1748", "1751", "1766", "1778", "1802",
"1815", "1817", "1819", "1828", "1830", "1831", "1834", "1835",
"1836", "1837", "1838", "1839", "1840", "1841", "1842", "1843",
"1844", "1845", "1846", "1847", "3097", "3164", "3168", "3169",
"3170", "3178", "3180", "3181", "3182"), class = "factor"), BLOCK = structure(c(12L,
77L, 163L, 67L, 28L, 170L, 90L, 36L, 52L, 2L, 15L, 19L, 168L,
103L, 188L, 31L, 203L, 66L, 29L, 46L, 34L, 32L, 27L, 16L, 83L,
48L, 82L, 30L, 171L, 14L, 115L, 54L, 93L, 65L, 50L, 187L, 58L,
91L, 200L, 6L, 169L, 135L, 99L, 148L, 101L, 104L, 107L, 128L,
153L, 146L, 41L, 22L, 53L, 87L, 131L, 151L, 110L, 10L, 44L, 11L,
13L, 20L, 42L, 202L, 111L, 38L, 183L, 51L, 199L, 109L, 75L, 134L,
92L, 166L, 182L, 97L, 100L, 1L, 86L, 181L, 25L, 108L, 94L, 116L,
72L, 18L, 23L, 76L, 185L, 81L, 62L, 63L, 56L, 204L, 85L, 95L,
129L, 49L, 147L, 106L, 145L, 205L, 73L, 207L, 105L, 24L, 43L,
8L, 167L, 164L, 3L, 96L, 184L, 45L, 74L, 39L, 89L, 4L, 152L,
130L, 165L, 40L, 57L, 70L, 206L, 186L, 7L, 37L, 9L, 102L, 132L,
127L, 88L, 80L, 98L, 139L, 196L, 174L, 118L, 215L, 194L, 193L,
208L, 172L, 122L, 143L, 141L, 123L, 161L, 209L, 213L, 178L, 159L,
160L, 191L, 177L, 192L, 144L, 175L, 211L, 140L, 180L, 173L, 125L,
119L, 120L, 210L, 214L, 136L, 154L, 162L, 190L, 158L, 216L, 142L,
124L, 212L, 195L, 155L, 121L, 64L, 68L, 117L, 59L, 71L, 35L,
69L, 201L, 21L, 84L, 61L, 114L, 17L, 112L, 55L, 150L, 113L, 79L,
78L, 47L, 33L, 149L, 60L, 189L, 5L, 133L, 26L, 137L, 197L, 179L,
126L, 198L, 157L, 176L, 138L, 156L), .Label = c("101", "102",
"103", "104", "105", "106", "107", "108", "109", "110", "111",
"112", "113", "114", "115", "116", "117", "118", "201", "202",
"203", "204", "205", "206", "207", "208", "209", "210", "211",
"212", "213", "214", "215", "216", "217", "218", "301", "302",
"303", "304", "305", "306", "307", "308", "309", "310", "311",
"312", "313", "314", "315", "316", "317", "318", "401", "402",
"403", "404", "405", "406", "407", "408", "409", "410", "411",
"412", "413", "414", "415", "416", "417", "418", "501", "502",
"503", "504", "505", "506", "507", "508", "509", "510", "511",
"512", "513", "514", "515", "516", "517", "518", "601", "602",
"603", "604", "605", "606", "607", "608", "609", "610", "611",
"612", "613", "614", "615", "616", "617", "618", "701", "702",
"703", "704", "705", "706", "707", "708", "709", "710", "711",
"712", "713", "714", "715", "716", "717", "718", "801", "802",
"803", "804", "805", "806", "807", "808", "809", "810", "811",
"812", "813", "814", "815", "816", "817", "818", "901", "902",
"903", "904", "905", "906", "907", "908", "909", "910", "911",
"912", "913", "914", "915", "916", "917", "918", "1001", "1002",
"1003", "1004", "1005", "1006", "1007", "1008", "1009", "1010",
"1011", "1012", "1013", "1014", "1015", "1016", "1017", "1018",
"1101", "1102", "1103", "1104", "1105", "1106", "1107", "1108",
"1109", "1110", "1111", "1112", "1113", "1114", "1115", "1116",
"1117", "1118", "1201", "1202", "1203", "1204", "1205", "1206",
"1207", "1208", "1209", "1210", "1211", "1212", "1213", "1214",
"1215", "1216", "1217", "1218"), class = "factor"), PLOT = structure(c(3L,
1L, 2L, 3L, 3L, 2L, 3L, 3L, 3L, 1L, 3L, 1L, 2L, 3L, 2L, 3L, 2L,
3L, 3L, 3L, 3L, 3L, 1L, 3L, 3L, 3L, 3L, 3L, 2L, 3L, 2L, 3L, 1L,
3L, 3L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 3L, 3L, 3L, 2L, 2L,
2L, 1L, 1L, 3L, 3L, 2L, 2L, 2L, 3L, 1L, 3L, 3L, 1L, 1L, 2L, 2L,
1L, 2L, 3L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 3L, 1L, 3L, 2L, 1L,
3L, 1L, 2L, 3L, 3L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 3L, 1L, 2L,
3L, 2L, 3L, 2L, 2L, 1L, 2L, 3L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L,
1L, 1L, 1L, 3L, 1L, 2L, 2L, 2L, 1L, 1L, 3L, 2L, 2L, 1L, 1L, 1L,
3L, 2L, 2L, 3L, 1L, 1L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 3L, 3L, 2L, 1L, 3L, 3L, 3L, 2L, 1L, 3L, 1L, 2L, 3L,
2L, 1L, 2L, 2L, 1L, 1L, 3L, 3L, 2L, 1L, 2L, 1L, 2L, 1L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("1", "2", "3", "4"), class = "factor"),
RANGE = structure(c(1L, 5L, 10L, 4L, 2L, 10L, 5L, 2L, 3L,
1L, 1L, 2L, 10L, 6L, 11L, 2L, 12L, 4L, 2L, 3L, 2L, 2L, 2L,
1L, 5L, 3L, 5L, 2L, 10L, 1L, 7L, 3L, 6L, 4L, 3L, 11L, 4L,
6L, 12L, 1L, 10L, 8L, 6L, 9L, 6L, 6L, 6L, 8L, 9L, 9L, 3L,
2L, 3L, 5L, 8L, 9L, 7L, 1L, 3L, 1L, 1L, 2L, 3L, 12L, 7L,
3L, 11L, 3L, 12L, 7L, 5L, 8L, 6L, 10L, 11L, 6L, 6L, 1L, 5L,
11L, 2L, 6L, 6L, 7L, 4L, 1L, 2L, 5L, 11L, 5L, 4L, 4L, 4L,
12L, 5L, 6L, 8L, 3L, 9L, 6L, 9L, 12L, 5L, 12L, 6L, 2L, 3L,
1L, 10L, 10L, 1L, 6L, 11L, 3L, 5L, 3L, 5L, 1L, 9L, 8L, 10L,
3L, 4L, 4L, 12L, 11L, 1L, 3L, 1L, 6L, 8L, 8L, 5L, 5L, 6L,
8L, 11L, 10L, 7L, 12L, 11L, 11L, 12L, 10L, 7L, 8L, 8L, 7L,
9L, 12L, 12L, 10L, 9L, 9L, 11L, 10L, 11L, 8L, 10L, 12L, 8L,
10L, 10L, 7L, 7L, 7L, 12L, 12L, 8L, 9L, 9L, 11L, 9L, 12L,
8L, 7L, 12L, 11L, 9L, 7L, 4L, 4L, 7L, 4L, 4L, 2L, 4L, 12L,
2L, 5L, 4L, 7L, 1L, 7L, 4L, 9L, 7L, 5L, 5L, 3L, 2L, 9L, 4L,
11L, 1L, 8L, 2L, 8L, 11L, 10L, 7L, 11L, 9L, 10L, 8L, 9L), .Label = c("1",
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12"), class = "factor"),
ROW = structure(c(12L, 5L, 1L, 13L, 10L, 8L, 18L, 18L, 16L,
2L, 15L, 1L, 6L, 13L, 8L, 13L, 5L, 12L, 11L, 10L, 16L, 14L,
9L, 16L, 11L, 12L, 10L, 12L, 9L, 14L, 7L, 18L, 3L, 11L, 14L,
7L, 4L, 1L, 2L, 6L, 7L, 9L, 9L, 4L, 11L, 14L, 17L, 2L, 9L,
2L, 5L, 4L, 17L, 15L, 5L, 7L, 2L, 10L, 8L, 11L, 13L, 2L,
6L, 4L, 3L, 2L, 3L, 15L, 1L, 1L, 3L, 8L, 2L, 4L, 2L, 7L,
10L, 1L, 14L, 1L, 7L, 18L, 4L, 8L, 18L, 18L, 5L, 4L, 5L,
9L, 8L, 9L, 2L, 6L, 13L, 5L, 3L, 13L, 3L, 16L, 1L, 7L, 1L,
9L, 15L, 6L, 7L, 8L, 5L, 2L, 3L, 6L, 4L, 9L, 2L, 3L, 17L,
4L, 8L, 4L, 3L, 4L, 3L, 16L, 8L, 6L, 7L, 1L, 9L, 12L, 6L,
1L, 16L, 8L, 8L, 13L, 16L, 12L, 10L, 17L, 14L, 13L, 10L,
10L, 14L, 17L, 15L, 15L, 17L, 11L, 15L, 16L, 15L, 16L, 11L,
15L, 12L, 18L, 13L, 13L, 14L, 18L, 11L, 17L, 11L, 12L, 12L,
16L, 10L, 10L, 18L, 10L, 14L, 18L, 16L, 16L, 14L, 15L, 11L,
13L, 10L, 14L, 9L, 5L, 17L, 17L, 15L, 3L, 3L, 12L, 7L, 6L,
17L, 4L, 1L, 6L, 5L, 7L, 6L, 11L, 15L, 5L, 6L, 9L, 5L, 7L,
8L, 11L, 17L, 17L, 18L, 18L, 13L, 14L, 12L, 12L), .Label = c("1",
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12",
"13", "14", "15", "16", "17", "18"), class = "factor"), YIELD = c(7882.814724,
7641.976671, 7535.187491, 8462.821158, 6470.762695, 7086.39647,
7260.626003, 8374.363239, 8225.545799, 6870.562479, 7260.303179,
6472.786879, 6535.801894, 7335.468082, 8101.853381, 7544.810974,
5597.940891, 8144.903193, 8489.541356, 7420.247609, 8267.229308,
7388.809243, 8753.922873, 7675.2452, 7540.083649, 7459.719121,
7614.590404, 6910.577593, 7655.161236, 8086.00529, 6754.554032,
9141.060314, 7728.70075, 7210.881432, 8872.660416, 7341.942246,
8211.265337, 9030.218757, 8957.01212, 7134.079145, 8580.60533,
8901.807114, 9009.635596, 8972.04225, 8850.07798, 7244.08863,
9357.355395, 7693.962907, 9059.604638, 8115.135788, 8073.220877,
7694.865425, 7168.389384, 7931.776306, 8310.054831, 7743.358631,
7241.417998, 7887.710882, 8671.335868, 7900.074562, 7089.929401,
8252.964285, 8038.601576, 8749.99335, 7880.418003, 7227.593551,
9733.562528, 7715.095262, 6926.775409, 7770.203085, 9000.211927,
7808.710708, 8239.82626, 8252.964285, 9546.314331, 2801.654022,
7865.302917, 6472.037973, 11286.93314, 7698.702989, 8239.164252,
8391.871173, 7817.085477, 7987.7324, 8517.420004, 8286.027753,
8021.268999, 8605.836444, 8360.390812, 8408.648702, 6980.52271,
8484.391646, 7604.489488, 8047.32564, 6859.736888, 8211.744547,
8338.224508, 7549.875965, 7831.170315, 8002.372075, 8092.398475,
7233.303386, 7880.198456, 6431.676768, 8146.454012, 9012.217125,
7696.760712, 7916.314754, 8372.430545, 4552.305881, 4744.119616,
8072.706265, 8038.601576, 8070.612573, 7631.800415, 8124.412039,
7958.686488, 8565.578204, 7204.2532, 7782.851494, 8195.743097,
8075.444598, 7468.681342, 7376.4572, 7019.132415, 7450.186973,
7900.853201, 7077.396698, 6781.366002, 8195.304822, 7581.211378,
8155.600681, 7446.611537, 7887.710882, 6849.690117, 6384.206298,
6965.647058, 7732.576444, 7687.296996, 7887.710882, 8061.034883,
7861.831189, 6690.298381, 7982.777954, 8310.054831, 7476.530867,
5840.137517, 8012.816166, 9211.484507, 8906.076566, 7227.155276,
6795.608201, 6926.023806, 8026.998142, 7388.809243, 7700.812705,
7493.134187, 7397.470718, 6794.411986, 8475.249868, 8387.892097,
8503.435859, 7890.106874, 7631.800415, 8349.757061, 7852.912013,
7758.848165, 7580.919692, 6402.21648, 6920.804051, 8628.194894,
7489.137138, 7866.037678, 7311.596266, 8746.497033, 9147.374207,
9022.033508, 8475.348448, 8911.007949, 8961.95446, 8476.003123,
8932.837953, 8661.336305, 8949.625535, 9048.100379, 10684.87284,
8845.185424, 8182.999872, 8986.675848, 8136.137692, 10504.2443,
8848.254372, 7233.813327, 8707.732966, 8381.547529, 10471.33626,
7682.888263, 8071.666541, 7428.171461, 9736.360333, 9378.789551,
8294.552055, 8225.545799, 8874.930993, 8459.226077, 8749.99335,
9192.455984, 7875.820212, 8982.410256, 8642.199262, 8935.14394,
8480.821358, 10240.80452, 8746.68483, 7619.897735, 8417.475201
)), .Names = c("ENTRY", "BLOCK", "PLOT", "RANGE", "ROW",
"YIELD"), row.names = 372:587, class = "data.frame")
Пространственное расположение данных:
library(reshape2)
dcast(my.data, RANGE ~ ROW, value.var ="YIELD")
Возможные примеры моделей для анализа данных:
library(nlme)
fit1 = lme(fixed = YIELD ~ ENTRY, data = my.data,
random= ~1 | BLOCK,
method = "ML")
fit2 = lme(fixed = YIELD ~ ENTRY, data = my.data,
random= ~1 | BLOCK,
corr = corSpatial(form = ~RANGE+ROW),
method = "ML")