HTML

Biostatblog

Biostatisztika a gyakorlatban kicsiknek és nagyoknak. Érdeklődőknek valamint kekeckedő profiknak. A közös nyelv amin beszélünk - a csodálatos R.

Friss topikok

Linkblog

Archívum

Illik ez az illesztés?! (2)

2012.02.04. 11:52 serpico

...megvan az adattömb, lássuk az ornátust! Azaz ábrázoljuk az árut a lehető legbénább módon, viszont gyorsan.

Col <- with(LongSubData,ifelse(Type=="n","red","blue"))
Pch <- with(LongSubData,ifelse(time=="Before",21,19))
plot(Weight~Days,col=Col,pch=Pch,data=LongSubData,cex.axis=2,cex.lab=2,cex=2)
legend("topleft",legend=c("Pucéran","Ruhában etetés előtt","Ruhában etetés után"),
pch=c(21,21,19),col=c("red","blue","blue"),cex=2)

Post02_Fig001.jpg


Minden növekedési görbéről tudvalevő – ok, a legtöbbről - hogy valamilyen felső határértékhez tart. (Szépemlékű irodalom tanárnőm állítása szerint ez a jelenség tántorította el a matematikától – folyamatosan közelíti, de soha el nem éri…az meg hogy lehet, he?). Ebben az esetben ezt zárójelbe tenném és azt feltételezném, hogy bár ismerjük e jelenséget, az adatok ezen a ponton még nem teszik lehetővé a korrekt aszimptotikus illesztést. Értsd; az adathalmaz lehetővé tesz ugyan ilyen illesztést, de úgyse hinnénk el mert ezek matematikai modellek, nem gondolkodnak csak illeszkednek – ha kell ha nem… Nekünk kell gondolkodni. Milyen szomorú! Maradjunk egyenlőre csak a lineáris közelítéseknél, azok olyan biztonságosak.
Kiindulásként három szisztematikus hatást feltételezek – ehhez nem is kell sok fantázia. Az idő előrehaladtával a súly nyilván növekszik, ha ruhában van az növeli a súlyt ha eszik az is. A lineáris modell a következőképpen néz ki:
lm01 <- lm(Weight~Days+time+Type,data=LongSubData)
> summary(lm01)
Call:
lm(formula = Weight ~ Days + time + Type, data = LongSubData)

Residuals:
     Min       1Q   Median       3Q      Max
-195.760  -52.532   -3.697   47.737  185.315

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3975.5414    16.0212  248.14   <2e-16 ***
Days          29.9103     0.3309   90.39   <2e-16 ***
timeBefore  -101.7033     7.3169  -13.90   <2e-16 ***
Types        181.3852    13.2069   13.73   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 69.8 on 393 degrees of freedom
Multiple R-squared: 0.9566,     Adjusted R-squared: 0.9563
F-statistic:  2888 on 3 and 393 DF,  p-value: < 2.2e-16

Bika kis modell mi? Lássuk mit is mond!
A gyerek 3976 g-al született…majdnem stimmel 4116 g volt. Az illesztett időszakban naponta 30 g-ot gyarapodott. A ráadott ruha súlya átlag 181 g-ot nyomott (nem páncélinget adtunk rá), valamint az etetések során 101 ml –t evett (hát nem teljesen…). Ezek mind olyan számok amik majdnem hihetőek, az egyetlen gyanús darab az a 70g-os RSE. Ez ha minden ok lenne a modellel a mérési hibát kéne hogy jelentse. Nem vagyok 100%-ig elégedett a termékkel de a legrosszabb időszakaiban is maximum a harmada volt ennek a hiba.
Mi lehet a baj? Több ötletem is van.
A gyerek napi súlyváltozása elméletileg inkább így néz ki
Post02_Fig002.jpg

Vagyis van egy napon belüli ritmusa a testsúlyingadozásnak. Ez már önmagában is jelentős nem magyarázott varianciát jelent. Kár hogy a mérések mellé rendelt sorszámok csak durva időfelbontást tennének lehetővé, plusz most csak lineáris modellekkel operálunk.
Van ezenfelül egy másik valószínű jelenség, ez pedig a megevett mennyiség változása az időben. Ez egy picit a feje tetejére fordul itt, mert a modellben ez úgy jelentkezne, hogy más meredekséget rendelünk az idő faktor mellé etetett és nem etetett esetben. Így ni:



lm02 <- lm(Weight~Days*time+Type,data=LongSubData)
>summary(lm02)
Call:
lm(formula = Weight ~ Days * time + Type, data = LongSubData)

Residuals:
     Min       1Q   Median       3Q      Max
-191.014  -49.101   -4.072   48.986  187.967

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     3956.1589    17.6821 223.738  < 2e-16 ***
Days              30.8026     0.4837  63.684  < 2e-16 ***
timeBefore       -65.0842    16.2752  -3.999  7.6e-05 ***
Types            181.0582    13.1190  13.801  < 2e-16 ***
Days:timeBefore   -1.6579     0.6593  -2.515   0.0123 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 69.33 on 392 degrees of freedom
Multiple R-squared: 0.9573,     Adjusted R-squared: 0.9569
F-statistic:  2197 on 4 and 392 DF,  p-value: < 2.2e-16

Na, ettől jobb is lett a modell.
Ennek bizonyítására összehasonlítom a két modellt.
> anova(lm01,lm02)
Analysis of Variance Table

Model 1: Weight ~ Days + time + Type
Model 2: Weight ~ Days * time + Type
  Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
1    393 1914670                              
2    392 1884274  1     30396 6.3235 0.01231 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Szóval ez a modell pont annyival jobb, mint amennyit az interakció hozzátesz (de meglepő…).
Az új modell szerint 3956 grammal született. A ruha még mindig 181 gramm DE, eleinte 65 mili körül eszik és ez minden nappal 1,6 milivel nő. Nem mondom hogy tökéletes, de egynek jó lesz.
Lássuk hogy grafikusan milyen az illesztett modell.
Post02_Fig003.jpgCsudás...


Szólj hozzá!

Címkék: ábra curve lieáris modell lm() plot() summary() legend()

Illik ez az illesztés?! (1)

2012.02.01. 12:52 serpico

Mióta a lányom "kipottyant", a feleségemmel elhitettem hogy ÉLETBEVÁGÓan fontos a gyerek súlyát naponta mérni. Egy ideig még az etetések előtt és utáni mérésre is sikerült rávenni, amelyek természetesen ugyancsak életbevágóan fontosak voltak, mind a magunk mind lányunk boldogsága és jóléte szempontjából.

A háttérben természetesen csak a bennem lévő monoton hang dolgozott; ADAT...ADAT...gyűljön! gyűljön! A pompás tervnek kb. a második "stagnálós" időszak vetett véget.

-Ez a gyerek nem hajlandó enni! És kurvára idegesít hogy ezt még látom is (minden este update-eltem neki a kiváló ábrát rövid elemzésekkel tűzdelve a grafikon profilját...)

Szóval adva van itt ez az adattömb. A kialakításáról annyit hogy akkor ez jó ötletnek tűnt...

A Type változó jelöli hogy a gyerek ruhában vagy mezítelenül lett e mérve, az Obs az adott napon belüli sorrendre utal, a Resp1 és 2 pedig a vonatkozó etetés előtti (Resp1) vagy utáni (Resp2) súlyok x*10 grammban kifejezve.

Az eredeti adattömböt Excel-ből olvastam be az RODBC csomag segítségével. Itt jegyezném meg hogy a biostatisztikus állandó problémája az, hogy a kutatók a táblázatkezelők miatt hajlamosak a "rövid-széles" (a finom erotikára fogékonyak itt kicsit elmélázhatnak...) adattömb kialakításra. Ennek oka nem az a bornírt szándék hogy minket megszívassanak, hanem az hogy a táblázatkezelőkbe bele lett hekkelve az ábra készítés funkció is - ami praktikusan képtelen a normális, "hosszú - keskeny" adattömb kezelésére. Így aztán a legváltozatosabb adattömb kialakításokkal lehet találkozni. Személyes kedvenceim az ún. színezők (kutyatenyésztők most ámulnak egyet). E típus az Excel cellaszíneit is felhasználják adataik jó esetben szépítésére de inkább valamilyen nem is elhanyagolható okból történő megkülönböztetésére és lelkesen magyarázzák hogy ez milyen jó! Nekik. Nekem HORROR. Az ilyen adatok gatyába rázása egy minőségbiztosítási rémálom (erről majd talán máskor írok egy pár sort...).

Visszakanyarodva, így utólag annyira már nem jó ötlet, kezdve azzal hogy a meztelen súly esetében (Type=="n") éppen fordítva van a bevitel, mivel itt kvázi mindig etetés előtti súlyról van szó. Ráadásul, ez a formátum a súly szempontjából egy részlegesen "rövid-széles" adattömbnek minősül. Először tehát ezen kéne segíteni!

R-ben erre a kiváló reshape() függvény használatos (SASolók mindezt behelyettesíthetik egy PROC TRANSPOSE-val).

LongData <- reshape(BabyWeight,varying = list(4:5),idvar = c("Date","Type","Obs"),
v.names = "Weight",direction = "long",times=c("Before","After"))

ilyen "szép" lett:

                            Date Type Obs Weight Days   time
2011-10-10.s.1.Before 2011-10-10    s   1    411    4 Before
2011-10-10.s.2.Before 2011-10-10    s   2    429    4 Before
2011-10-10.s.3.Before 2011-10-10    s   3    409    4 Before
2011-10-10.s.4.Before 2011-10-10    s   4    428    4 Before

A fura izék (pl.: 2011-10-10.s.1.Before) az első "oszlopban" a sornevek.

Akkor megvan az adattömb. Hurrá! Kár hogy van egy rahedli nullás rekord, néhány hibás adatbevitel meg hogy az etetés előtti-utáni méricskézés valahol a 40. naptól almás...szóval szűrünk.

LongSubData <- subset(LongData, Days<40 & Weight>100)

Aztán kijavítjuk a "pucér" mérés típus hibáját
LongSubData$time[with(LongSubData,Type=="n")] <- "Before"

Valamint a babamérleg 10g-os pontosságig jelez ki tehát:
LongSubData$Weight <- with(LongSubData,Weight*10)

Na ennyit az eleinte jó ötletnek tűnt típusú adattömbről - dióhéjban.

Szólj hozzá!

Címkék: subset() RODBC reshape()

süti beállítások módosítása