Próbuję stworzyć postać, która pokazuje związek między kopiami wirusów a pokryciem genomu (GCC). Tak wyglądają moje dane:
Na początku po prostu nakreśliłem regresję liniową, ale moi przełożeni powiedzieli mi, że to nieprawda, i wypróbowałem krzywą sigmoidalną. Zrobiłem to za pomocą geom_smooth:
library(scales)
ggplot(scatter_plot_new, aes(x = Copies_per_uL, y = Genome_cov, colour = Virus)) +
geom_point() +
scale_x_continuous(trans = log10_trans(), breaks = trans_breaks("log10", function(x) 10^x), labels = trans_format("log10", math_format(10^.x))) +
geom_smooth(method = "gam", formula = y ~ s(x), se = FALSE, size = 1) +
theme_bw() +
theme(legend.position = 'top', legend.text = element_text(size = 10), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size=12), axis.title.y = element_text(margin = margin (r = 10)), axis.title.x = element_text(margin = margin(t = 10))) +
labs(x = "Virus copies/µL", y = "GCC (%)") +
scale_y_continuous(breaks=c(25,50,75,100))
Jednak moi przełożeni twierdzą, że jest to również nieprawidłowe, ponieważ krzywe sprawiają, że wygląda na to, że GCC może przekroczyć 100%, czego nie może.
Moje pytanie brzmi: jaki jest najlepszy sposób na pokazanie związku między kopiami wirusów a GCC? Chcę wyjaśnić, że A) niska liczba kopii wirusa = niski GCC i że B) po określonej ilości wirusa kopiuje płaskowyże GCC.
Badałem wiele różnych metod - GAM, LOESS, logistyka, fragmentarycznie - ale nie wiem, jak powiedzieć, która metoda jest najlepsza dla moich danych.
EDYCJA: to są dane:
>print(scatter_plot_new)
Subsample Virus Genome_cov Copies_per_uL
1 S1.1_RRAV RRAV 100 92500
2 S1.2_RRAV RRAV 100 95900
3 S1.3_RRAV RRAV 100 92900
4 S2.1_RRAV RRAV 100 4049.54
5 S2.2_RRAV RRAV 96.9935 3809
6 S2.3_RRAV RRAV 94.5054 3695.06
7 S3.1_RRAV RRAV 3.7235 86.37
8 S3.2_RRAV RRAV 11.8186 84.2
9 S3.3_RRAV RRAV 11.0929 95.2
10 S4.1_RRAV RRAV 0 2.12
11 S4.2_RRAV RRAV 5.0799 2.71
12 S4.3_RRAV RRAV 0 2.39
13 S5.1_RRAV RRAV 4.9503 0.16
14 S5.2_RRAV RRAV 0 0.08
15 S5.3_RRAV RRAV 4.4147 0.08
16 S1.1_UMAV UMAV 5.7666 1.38
17 S1.2_UMAV UMAV 26.0379 1.72
18 S1.3_UMAV UMAV 7.4128 2.52
19 S2.1_UMAV UMAV 21.172 31.06
20 S2.2_UMAV UMAV 16.1663 29.87
21 S2.3_UMAV UMAV 9.121 32.82
22 S3.1_UMAV UMAV 92.903 627.24
23 S3.2_UMAV UMAV 83.0314 615.36
24 S3.3_UMAV UMAV 90.3458 632.67
25 S4.1_UMAV UMAV 98.6696 11180
26 S4.2_UMAV UMAV 98.8405 12720
27 S4.3_UMAV UMAV 98.7939 8680
28 S5.1_UMAV UMAV 98.6489 318200
29 S5.2_UMAV UMAV 99.1303 346100
30 S5.3_UMAV UMAV 98.8767 345100
źródło
method.args=list(family=quasibinomial))
argumenty dogeom_smooth()
oryginalnego kodu ggplot.se=FALSE
. Zawsze miło jest pokazywać ludziom, jak duża jest niepewność ...Odpowiedzi:
Innym sposobem na obejście tego byłoby użycie sformułowania bayesowskiego, może być nieco ciężkie na początku, ale zwykle znacznie ułatwia wyrażanie specyfiki twojego problemu, a także uzyskiwanie lepszych pomysłów na to, gdzie „niepewność” jest
Stan jest samplerem Monte Carlo ze stosunkowo łatwym w użyciu interfejsem programowym, biblioteki są dostępne dla R i innych, ale tutaj używam Pythona
używamy sigmoidu jak wszyscy: ma motywacje biochemiczne, a także jest matematycznie bardzo wygodny w pracy. niezłą parametryzacją tego zadania jest:
gdzie
alpha
określa punkt środkowy krzywej sigmoidalnej (tj. gdzie przecina 50%) ibeta
określa nachylenie, wartości bliższe zeru są bardziej płaskieaby pokazać, jak to wygląda, możemy pobrać twoje dane i wykreślić je za pomocą:
gdzie
raw_data.txt
zawiera dane, które podałeś, a ja przekształciłem zasięg w coś bardziej użytecznego. współczynniki 5.5 i 3 wyglądają ładnie i dają wykres bardzo podobny do innych odpowiedzi:aby „dopasować” tę funkcję za pomocą Stana, musimy zdefiniować nasz model za pomocą własnego języka będącego mieszanką R i C ++. prosty model byłby mniej więcej taki:
co mam nadzieję, że brzmi OK. mamy
data
blok, który definiuje dane, których oczekujemy podczas próbkowania modelu,parameters
definiuje rzeczy, które są próbkowane, imodel
definiuje funkcję prawdopodobieństwa. Mówisz Stanowi, aby „skompilował” model, co zajmuje chwilę, a następnie możesz próbkować z niego przy użyciu niektórych danych. na przykład:arviz
ułatwia tworzenie ładnych wykresów diagnostycznych, a drukowanie pasowania daje ładne podsumowanie parametrów w stylu R:duże odchylenie standardowe
beta
mówi, że dane naprawdę nie dostarczają zbyt wielu informacji na temat tego parametru. również niektóre z odpowiedzi, które wskazują ponad 10 cyfr znaczących w modelu, są nieco zawyżoneponieważ w niektórych odpowiedziach zauważono, że każdy wirus może potrzebować własnych parametrów, rozszerzyłem model, aby zezwolić
alpha
ibeta
różnić się w zależności od „wirusa”. wszystko robi się trochę dziwnie, ale dwa wirusy prawie na pewno mają różnealpha
wartości (tj. potrzebujesz więcej kopii / μL RRAV dla tego samego zasięgu), a wykres pokazujący to:dane są takie same jak poprzednio, ale narysowałem krzywą dla 40 próbek tylnej.
UMAV
wydaje się stosunkowo dobrze określony, chociażRRAV
może podążać tym samym nachyleniem i wymagać większej liczby kopii lub mieć bardziej strome nachylenie i podobną liczbę kopii. większa część tylnej masy wymaga większej liczby kopii, ale ta niepewność może wyjaśnić niektóre różnice w innych odpowiedziach dotyczących znalezienia różnych rzeczyNajczęściej stosowane odpowiadając to jako ćwiczenie poprawić moją znajomość Stan, a ja umieścić Jupyter zeszyt to tu w przypadku gdy ktoś jest zainteresowany / chce powtórzyć to.
źródło
(Edytowane z uwzględnieniem komentarzy poniżej. Podziękowania dla @BenBolker i @WeiwenNg za pomocne informacje).
Dopasuj ułamkową regresję logistyczną do danych. Jest dobrze dostosowany do danych procentowych, które są ograniczone od 0 do 100% i jest dobrze uzasadniony teoretycznie w wielu obszarach biologii.
Należy pamiętać, że może być konieczne podzielenie wszystkich wartości przez 100, aby je dopasować, ponieważ programy często oczekują, że dane będą zawierać się w przedziale od 0 do 1. I, jak zaleca Ben Bolker, aby rozwiązać ewentualne problemy spowodowane ścisłymi założeniami rozkładu dwumianowego dotyczącymi wariancji, użyj zamiast tego rozkład quasibinomial.
Podjąłem pewne założenia na podstawie twojego kodu, na przykład, że interesują Cię 2 wirusy, które mogą wykazywać różne wzorce (tj. Może istnieć interakcja między typem wirusa a liczbą kopii).
Po pierwsze, model pasuje:
Jeśli ufasz wartościom p, dane wyjściowe nie sugerują, że oba wirusy różnią się znacząco. Jest to sprzeczne z poniższymi wynikami @ NickCox, chociaż zastosowaliśmy różne metody. W każdym razie nie byłbym bardzo pewny siebie z 30 punktami danych.
Po drugie, spisek:
Nie jest trudno samodzielnie kodować sposób wizualizacji wyników, ale wydaje się, że istnieje pakiet ggPredict, który wykona większość pracy za ciebie (nie mogę za to ręczyć, sam tego nie próbowałem). Kod będzie wyglądał mniej więcej tak:Aktualizacja: Nie polecam już bardziej ogólnie kodu ani funkcji ggPredict. Po wypróbowaniu okazało się, że wykreślone punkty nie odzwierciedlają dokładnie danych wejściowych, ale zostały zmienione z jakiegoś dziwnego powodu (niektóre z wykreślonych punktów były powyżej 1 i poniżej 0). Zalecam więc samodzielne kodowanie, ale to więcej pracy.
źródło
family=quasibinomial()
aby uniknąć ostrzeżenia (i podstawowych problemów związanych ze zbyt ścisłymi założeniami dotyczącymi wariancji). Skorzystaj z porady @ mkt na temat drugiego problemu.To nie jest inna odpowiedź niż @mkt, ale w szczególności wykresy nie mieszczą się w komentarzu. Najpierw dopasowuję krzywą logistyczną w Stata (po zalogowaniu predyktora) do wszystkich danych i otrzymuję ten wykres
Równanie jest
100
invlogit
(-4,192654 + 1,880951log10
(Copies
))Teraz dopasowuję krzywe osobno dla każdego wirusa w najprostszym scenariuszu definiowania zmiennej wskaźnikowej przez wirusa. Tutaj do nagrania jest skrypt Stata:
To mocno naciska na niewielki zestaw danych, ale wartość P dla wirusa wydaje się wspierać dopasowanie dwóch krzywych razem.
źródło
Wypróbuj funkcję sigmoidalną . Istnieje wiele formuł tego kształtu, w tym krzywa logistyczna. Styczna hiperboliczna to kolejny popularny wybór.
Biorąc pod uwagę wykresy, nie mogę również wykluczyć prostej funkcji kroku. Obawiam się, że nie będziesz w stanie odróżnić funkcji kroku od dowolnej liczby specyfikacji sigmoid. Nie masz żadnych obserwacji, w których procent mieści się w zakresie 50%, więc prosty krok może być najbardziej oszczędnym wyborem, który nie ustępuje bardziej niż bardziej złożone modele
źródło
Oto dopasowania 4PL (logistyka 4 parametrów), zarówno ograniczone, jak i nieograniczone, z równaniem według CA Holsteina, M. Griffina, J. Honga, PD Sampsona, „Statystyczna metoda określania i porównywania granic wykrywania testów biologicznych”, Anal . Chem. 87 (2015) 9795–9801. Równanie 4PL pokazano na obu rysunkach, a znaczenie parametru jest następujące: a = niższa asymptota, b = współczynnik nachylenia, c = punkt przegięcia, d = górna asymptota.
Rycina 1 ogranicza a do równości 0% id do równości 100%:
Rycina 2 nie ma ograniczeń dotyczących 4 parametrów w równaniu 4PL:
To była świetna zabawa, nie udaję, że wiem coś biologicznego i ciekawe będzie, jak to wszystko się ułoży!
źródło
Wyodrębniłem dane z twojego wykresu rozrzutu, a moje wyszukiwanie równań ujawniło 3-parametrowe równanie typu logistycznego jako dobrego kandydata: „y = a / (1.0 + b * exp (-1.0 * c * x))”, gdzie „ x ”to podstawa dziennika 10 na twoją działkę. Dopasowane parametry to = = 9.0005947126706630E + 01, b = 1,2831794858584102E + 07, oraz c = 6,6483431489473155E + 00 dla moich wyodrębnionych danych, dopasowanie oryginalnych danych (log 10 x) powinno dać podobne wyniki, jeśli ponownie dopasujesz oryginalne dane, używając moich wartości jako wstępnych oszacowań parametrów. Moje wartości parametrów dają R-kwadrat = 0,983 i RMSE = 5,625 na wyodrębnionych danych.
EDYCJA: Teraz, gdy pytanie zostało poddane edycji w celu uwzględnienia rzeczywistych danych, oto wykres wykorzystujący powyższe równanie 3-parametrowe i szacunkowe parametry początkowe.
źródło
Ponieważ musiałem otworzyć swoje wielkie usta na temat Heaviside, oto wyniki. Ustawiłem punkt przejścia na log10 (viruscopies) = 2,5. Następnie obliczyłem standardowe odchylenia dwóch połówek zbioru danych - to znaczy, że Heaviside zakłada, że dane po obu stronach mają wszystkie pochodne = 0.
Od lewej strony std dev = 4,76 Po
lewej stronie std dev = 7,72
Ponieważ okazuje się, że w każdej partii jest 15 próbek, ogólny std dev jest średnią, czyli 6,24.
Zakładając, że „RMSE” cytowany w innych odpowiedziach to ogólnie „błąd RMS”, wydaje się, że funkcja Heaviside działa co najmniej tak dobrze, a jeśli nie lepiej niż, większość „krzywej Z” (zapożyczonej z nomenklatury odpowiedzi fotograficznej) pasuje tutaj.
edytować
Bezużyteczny wykres, ale wymagany w komentarzach:
źródło