Chcę dopasować model mieszany za pomocą pakietu lme4, nme, pakietu regresji baysian lub dowolnego dostępnego.
Model mieszany w konwencjach kodowania Asreml-R
zanim przejdziemy do szczegółów, możemy chcieć uzyskać szczegółowe informacje na temat konwencji asreml-R dla tych, którzy nie znają kodów ASREML.
y = Xτ + Zu + e ........................(1) ; zwykły model mieszany z, y oznacza wektor obserwacji n × 1, gdzie τ jest wektorem ustalonych efektów p × 1, X jest macierzą projektową n × p o pełnej randze kolumny, która łączy obserwacje z odpowiednią kombinacją ustalonych efektów , u jest wektorem losowych efektów q × 1, Z jest macierzą projektową n × q, która łączy obserwacje z odpowiednią kombinacją efektów losowych, a e jest wektorem błędów resztkowych n × 1. Model (1) nazywa się liniowy model mieszany lub liniowy model mieszanych efektów. Zakłada się
gdzie macierze G i R są funkcjami odpowiednio parametrów γ i φ.
Parametr θ jest parametrem wariancji, który będziemy określać jako parametr skali.
W mieszanych modelach efektów z więcej niż jedną wariancją rezydualną, powstającą na przykład w analizie danych z więcej niż jedną sekcją lub zmienną, parametr θ jest ustawiony na jeden. W mieszanych modelach efektów z pojedynczą wariancją rezydualną wtedy θ jest równa wariancji rezydualnej (σ2). W takim przypadku R musi być macierzą korelacji. Dalsze szczegóły dotyczące modeli znajdują się w podręczniku Asreml (link) .
Struktury wariancji dla błędów: struktura R i struktury wariancji dla efektów losowych: można określić struktury G.
modelowanie wariancji w asreml () ważne jest zrozumienie tworzenia struktur wariancji za pomocą produktów bezpośrednich. Zwykle założeniem najmniejszych kwadratów (i domyślnym w asreml ()) jest to, że są one niezależnie i identycznie rozmieszczone (IID). Jeśli jednak dane pochodzą z eksperymentu polowego ułożonego w prostokątną tablicę r wierszy po c kolumnach, powiedzmy, moglibyśmy uporządkować reszty e jako macierz i potencjalnie uznać, że były one autokorelowane w wierszach i kolumnach. wektor w kolejności pól, tzn. sortując rzędy reszt w kolumnach (wykresy w blokach), wariancja reszt może być wtedy
są odpowiednio macierzami korelacji dla modelu wiersza (rząd r, parametr autokorelacji ½r) i modelu kolumny (rząd c, parametr autokorelacji ½c). Mówiąc dokładniej, dla typowych błędów w analizie prób terenowych zakłada się czasami dwuwymiarową, separowalną autoregresyjną strukturę przestrzenną (AR1 x AR1).
Przykładowe dane:
nin89 pochodzi z biblioteki asreml-R, gdzie hodowano różne odmiany w replikacjach / blokach w polu prostokątnym. Aby kontrolować dodatkową zmienność w kierunku wiersza lub kolumny, każdy wykres jest określany jako zmienne wiersza i kolumny (projekt kolumny wiersza). Tak więc ten projekt kolumny wiersza z blokowaniem. Wydajność jest mierzona zmiennie.
Przykładowe modele
Potrzebuję czegoś równoważnego kodom asreml-R:
Prosta składnia modelu będzie wyglądać następująco:
rcb.asr <- asreml(yield ∼ Variety, random = ∼ Replicate, data = nin89)
.....model 0
Model liniowy jest określony w argumentach stałych (wymaganych), losowych (opcjonalnych) i rcov (komponent błędu) jako obiekty formuły. Domyślnie jest to prosty termin błędu i nie trzeba go formalnie określać dla terminu błędu, jak w modelu 0 .
tutaj odmiana ma ustalony efekt, a losowa jest replikowana (bloki). Oprócz losowych i stałych terminów możemy określić termin błędu. Który jest domyślny w tym modelu 0. Resztkowy lub błąd komponentu modelu jest określony w obiekcie formuły za pomocą argumentu rcov, patrz następujące modele 1: 4.
Poniższy model1 jest bardziej złożony, w którym określono zarówno strukturę G (losową), jak i R (błąd).
Model 1:
data(nin89)
# Model 1: RCB analysis with G and R structure
rcb.asr <- asreml(yield ~ Variety, random = ~ idv(Replicate),
rcov = ~ idv(units), data = nin89)
Ten model jest równoważny powyższemu modelowi 0 i wprowadza użycie modelu wariancji G i R. W tym przypadku opcja random i rcov określa formuły random i rcov, aby jawnie określić struktury G i R. gdzie idv () to specjalna funkcja modelu w asreml (), która identyfikuje model wariancji. Wyrażenie idv (jednostki) jawnie ustawia macierz wariancji dla e na skalowaną tożsamość.
# Model 2: dwuwymiarowy model przestrzenny z korelacją w jednym kierunku
sp.asr <- asreml(yield ~ Variety, rcov = ~ Column:ar1(Row), data = nin89)jednostki eksperymentalne z nin89 są indeksowane według kolumn i rzędów. Oczekujemy więc losowej zmienności w dwóch kierunkach - w tym przypadku w rzędzie i kolumnie. gdzie ar1 () to specjalna funkcja określająca model wariancji autoregresji wariantu pierwszego rzędu dla Row. To wywołanie określa dwuwymiarową strukturę przestrzenną dla błędu, ale z korelacją przestrzenną tylko w kierunku wiersza. Model wariancji dla Kolumny to tożsamość (id ()), ale nie musi być formalnie określony, ponieważ jest to ustawienie domyślne.
# model 3: dwuwymiarowy model przestrzenny, struktura błędów w obu kierunkach
sp.asr <- asreml(yield ~ Variety, rcov = ~ ar1(Column):ar1(Row),
data = nin89)
sp.asr <- asreml(yield ~ Variety, random = ~ units,
rcov = ~ ar1(Column):ar1(Row), data = nin89)
podobny do powyższego modelu 2, jednak korelacja jest dwukierunkowa - autoregresyjna.
Nie jestem pewien, ile z tych modeli jest możliwych dzięki pakietom R. Nawet jeśli rozwiązanie któregokolwiek z tych modeli będzie bardzo pomocne. Nawet jeśli pakiet +50 może stymulować do opracowania takiego pakietu, będzie bardzo pomocny!
Zobacz MAYSaseen dostarczył dane wyjściowe z każdego modelu i dane (jako odpowiedź) do porównania.
Edycje: Oto sugestia, którą otrzymałem na forum dyskusyjnym modeli mieszanych: „Możesz spojrzeć na pakiety regresji i przestrzennej kowariancji Davida Clifforda. Ten pierwszy umożliwia dopasowanie modeli mieszanych (Gaussa), w których możesz bardzo elastycznie określać strukturę macierzy kowariancji (na przykład użyłem go do danych rodowodu). Pakiet spatialCovariance używa regresji, aby dostarczyć bardziej rozbudowanych modeli niż AR1xAR1, ale może mieć zastosowanie. Być może będziesz musiał porozmawiać z autorem o zastosowaniu go do konkretnego problemu. ”
lme4. Czy możesz (a) powiedzieć nam, dlaczego musisz to zrobić,lme4zamiastasreml-R(b) rozważyć opublikowanie w miejscu,r-sig-mixed-modelsgdzie jest bardziej odpowiednia wiedza specjalistyczna?corStructwnlme(dla korelacji anizotropowych) ... Byłoby pomocne, gdybyś mógł krótko określić (słowami lub równaniami) modele statystyczne odpowiadające tym stwierdzeniom ASREML, ponieważ nie wszyscy jesteśmy zaznajomieni z Składnia ASREML ...MCMCglmmi jestem tego pewien (pozaspatialCovariancewspomniano, którym jestem zaznajomiony z), że jedynym sposobem, aby to zrobić w R jest poprzez definiowanie nowychcorStructs - co jest możliwe, ale nie banalne.Odpowiedzi:
Możesz dopasować ten model do AD Model Builder. AD Model Builder to darmowe oprogramowanie do budowania ogólnych modeli nieliniowych, w tym ogólnych nieliniowych modeli efektów losowych. Na przykład można dopasować ujemny dwumianowy model przestrzenny, w którym zarówno średnia, jak i nadmierna dyspersja miały strukturę ar (1) x ar (1). Zbudowałem kod dla tego przykładu i dopasowałem go do danych. Jeśli ktoś jest zainteresowany, prawdopodobnie lepiej omówić to na liście na stronie http://admb-project.org
Uwaga: Istnieje wersja R ADMB, ale funkcje dostępne w pakiecie R są podzbiorem samodzielnego oprogramowania ADMB.
W tym przykładzie łatwiej jest utworzyć plik ASCII z danymi, wczytać go do programu ADMB, uruchomić program, a następnie wczytać oszacowania parametrów itp. Z powrotem do R dla tego, co chcesz zrobić.
Powinieneś zrozumieć, że ADMB nie jest zbiorem pakietów, ale raczej językiem do pisania oprogramowania do nieliniowej oceny parametrów. Jak powiedziałem wcześniej, lepiej omówić to na liście ADMB, gdzie wszyscy wiedzą o oprogramowaniu. Po zakończeniu i zrozumieniu modelu możesz opublikować wyniki tutaj. Oto link do kodów ML i REML, które zestawiłem dla danych dotyczących pszenicy.
http://lists.admb-project.org/pipermail/users/attachments/20111124/448923c8/attachment.zip
źródło
Model 0
ASReml-R
lme4
nlme
źródło
Model 1
ASReml-R
nlme
Zobacz sztuczkę
źródło
Model 2
ASReml-R
nlme
Pracuję, ale jeszcze nie rozgryzłem. Może być coś takiego. Nadal nie mógł dowiedzieć się, jak to zrobić
rcov=~Column:ar1(Row)znlmeźródło
Model 3
ASReml-R
nlme
Pracuję, ale jeszcze nie rozgryzłem. Może być coś takiego. Nadal nie mógł dowiedzieć się, jak to zrobić
rcov=~ar1(Column):ar1(Row)znlmeNie mogłem wymyślić, jak dopasować Model 2 i 3
nlme. Pracuję nad tym i zaktualizuję odpowiedź, kiedy to zrobisz. AleASReml-Rdo celów porównawczych podałem dane wyjściowe dla Modelu 2 i 3. Kevin ma duże doświadczenie w analizowaniu takich modeli, a Ben Bolker ma wspaniały autorytet w zakresie modeli mieszanych. Mam nadzieję, że mogą nam pomóc w modelach 2 i 3.źródło