Xử lý dữ liệu thống kê nông nghiệp với phần mềm excel và sas

Chia sẻ: bluesky_12

Cuốn tài liệu Xử lý dữ liệu thống kê nông nghiệp với phần mềm Excel và SAS được biên soạn nhằm phân tích dữ liệu từ các mô hình thiết kế thí thường gặp trong nghiên cứu nông nghiệp bằng phần mềm MS Excel và SAS. Tài liệu Xử lý dữ liệu thống kê nông nghiệp với phần mềm Excel và SAS được viết trong khuôn khổ Dự án công nghệ thông tin, Trường Đại học Nông nghiệp Hà Nội.

Bạn đang xem 10 trang mẫu tài liệu này, vui lòng download file gốc để xem toàn bộ.

Nội dung Text: Xử lý dữ liệu thống kê nông nghiệp với phần mềm excel và sas

 

  1. TRƯ NG ð I H C NÔNG NGHI P HÀ N I PGS. TS. Nguy n H i Thanh − ThS. ð ð c L c X LÝ D LI U TH NG KÊ NÔNG NGHI P v i ph n m m Excel và SAS (Bài gi ng cho d án CNTT) HÀ N I, THÁNG 10 NĂM 2008
  2. Ph n 1 X lý d li u th ng kê nông nghi p trong Excel Ph n m m Excel cho phép phân tích d li u nói chung, d li u sinh h c và nông nghi p nói riêng, m t cách khá hi u qu thông qua vi c s d ng menu Tools> Data Analysis (n u không có m c này thì ch n Tools> Add-in > Analysis ToolPak ñ cài ñ t thêm). Sau ñây là m t s công c x lý s li u th ng kê mà Excel cung c p. 1. Gi i thi u v phương pháp m u và th ng kê mô t 1.1. Gi i thi u v phương pháp kh o sát m u ð nghiên c u v m t ch s nào ñó trên các cá th c a m t t ng th v i r t nhi u cá th , có th ti n hành theo hai cách. Cách 1: ði u tra ch s ñó trên t t c các cá th c a t ng th ; cách này ñòi h i chi phí cao, t n kém th i gian, nhi u khi không kh thi. Cách 2: ði u tra m t m u ng u nhiên các cá th c a t ng th ; d a trên k t qu c a m u ñi u tra ñư c và các ñ nh lý c a lý thuy t xác su t c n x lý s li u m u ñ ñưa ra m t suy ñoán th ng kê v ch s ñó cho toàn b t ng th . Cách 2 có th ñư c g i là phương pháp th ng kê toán h c, hay còn g i là phương pháp kh o sát m u. T ng quát hơn, phương pháp kh o sát m u ñư c áp d ng khi c n nghiên c u m t s ch s nào ñó cũng như các m i liên quan c a chúng trên các cá th c a t ng th . Ký hi u X là ch s ng u nhiên mà chúng ta c n kh o sát trên các cá th c a m t t ng th . Xét m t m u ng u nhiên dung lư ng n c a X là (X1, X2, ..., Xn) trong ñó Xi, i = 1, 2, …, n, là các bi n ng u nhiên ñ c l p thu ñư c t X. ð ñơn gi n chúng ta g i m u này là m u lý thuy t. Tương ng v i m u lý thuy t trên là m u th c nghi m (x1, x2, ..., xn) trong ñó xi là giá tr ño ñư c c a Xi thu ñư c t k t qu c a th c nghi m. 1 ð i lư ng th ng kê X = (X1+X2 + ... + Xn) ñư c g i là trung bình m u lý thuy t n và ñư c l y làm ư c lư ng cho kỳ v ng E(X) c a X, E(X) ñư c coi là s ñ c trưng cho trung bình chung c a ch s X. ðây là m t ư c lư ng r t t t v i các tính ch t: không 1 ch ch, v ng và hi u qu . Lúc ñó, x = (x1+ x2+ ...+ xn) ñư c g i là trung bình m u th c n nghi m, chính là giá tr trung bình c a ch s X trên m u th c nghi m. Nh có tính ch t v ng c a ư c lư ng, khi dung lư ng m u khá l n, ñ l ch gi a trung bình chung và trung bình th c nghi m là khá nh trong h u h t các l n ti n hành th c nghi m. n ˆ1 ∑ (X − X) 2 ñư c g i là phương sai m u lý Tương t , ñ i lư ng th ng kê: S2 = i n i =1 1n ∑ (X i − X) 2 g i là phương sai thuy t chưa hi u ch nh, còn ñ i lư ng th ng kê S2 = n − 1 i =1 m u ñã hi u ch nh. Chúng ñ u ñư c l y làm ư c lư ng cho phương sai V(X) c a X v i V(X) ñư c coi là s ñ c trưng cho ñ bi n ñ ng c a ch s X xung quanh E(X). Trong khi 2
  3. phương sai m u lý thuy t chưa hi u ch nh ch có tính ch t v ng, thì phương sai m u lý thuy t ñã hi u ch nh l i có c ba tính ch t không ch ch, v ng và hi u qu . n 1 ∑ (x − x ) 2 cũng ñư c g i là phương sai m u th c nghi m chưa Lúc ñó, s 2 = ˆ i n i =1 n 1 ∑ (x − x ) 2 cũng ñư c g i là phương sai m u th c nghi m ñã hi u ch nh, còn s2 = n −1 i i =1 ˆ ˆ hi u ch nh tương ng v i m u th c nghiêm ñã có. S và s ñư c g i là các ñ l ch chu n m u lý thuy t và th c nghi m chưa hi u ch nh, còn S và s g i là ñ l ch chu n m u lý thuy t và th c nghi m ñã hi u ch nh. 1.2. Th ng kê mô t Sau khi có s li u m u th c nghi m, có th s d ng ch c năng th ng kê mô t trong phân tích s li u c a Excel ñ tính các s ñ c trưng m u c a m u th c nghi m như trung bình, ñ l ch chu n, sai s chu n, trung v , mode... S li u tính toán ñư c b trí theo c t ho c theo hàng. a. Các bư c th c hi n trong Excel Ch n Tools>Data Analysis>Descriptive Statistics, và khai báo các m c sau trong h p tho i: - Input range: mi n d li u k c nhãn. - Grouped by: Column (s li u theo c t). - Labels in first row : ðánh d u √ vào ô này n u có nhãn hàng ñ u. - Confidence level for mean: 95% ( ñ tin c y 95%). - K-th largest: 1 (1 S l n nh t, 2 s l n nhì ). - K-th smallest: 1 (1 S nh nh t, 2 s nh nhì ). - Output range: mi n ra. - Summary Statistics: ðánh d u √ n u mu n hi n các th ng kê cơ b n. Ví d 1: Kh o sát v các ñ c tính c a lúa ta thu ñư c b n c t s li u: dài bông (cm), P1000 (tr ng lư ng 1000 h t), s bông/m t cây, năng su t. S li u ñư c b trí như trong hình I.1. Ch n Tools>Data Analysis>Descriptive Statistics, sau ñó khai báo h p tho i. Hình I.1. B ng s li u kh o sát v lúa và khai báo h p tho i. 3
  4. K t qu thu ñư c cho trong hình I.2 Hình I.2. K t qu th ng kê mô t v các ñ c tính c a lúa b. Phân tích các k t qu thu ñư c M t s nh n xét sơ b trên các th ng kê thu ñư c như sau: - Mean cho ta giá tr trung bình c a dãy s . - Standard error cho ta bi t t s ñ l ch chu n m u /căn b c hai c a n. - Median cho giá tr ñi m gi a c a dãy s . N u 2 giá tr Mean và Median x p x nhau ta thì s li u là cân ñ i. Trong ví d 1 các c t s li u là cân ñ i tr c t “s bông” hơi b l ch. - Mode cho bi t giá tr x y ra nhi u nh t trên m u - Phương sai m u hay ñ l ch chu n m u (ñã hi u ch nh) cho ta bi t ñ phân tán c a s li u quanh giá tr trung bình, n u các giá tr này càng nh ch ng t s li u càng t p trung. - Kurtosis ñánh giá ñư ng m t ñ phân ph i c a dãy s li u có nh n hơn hay tù hơn ñư ng m t ñ chu n t c. N u trong kho ng t -2 ñ n 2 thì có th coi s li u x p x chu n. - Skewness ñánh giá ñư ng phân ph i l ch trái hay l ch ph i. N u trong kho ng t -2 ñ n 2 thì có th coi s li u cân ñ i g n như s li u trong phân ph i chu n t c. - Confidence Level ñư c hi u là n a ñ dài kho ng tin c y. Gi s Confidence Level là m thì kho ng tin c y c a trung bình t ng th là: (Mean- m , Mean+m). Trong ví d 1, hình I.2., ta có kho ng tin c y 95% c a “dài bông” là: ( 26.4- 0.55 , 26.4 +0.55), t c là (25.85 , 26.95). Trong trư ng h p m u có dung lư ng n không l n l m và phươ ng sai lý thuy t σ2 = V(X) chưa bi t, ta có công th c tìm kho ng tin c y v i ñ tin c y p = 1−α s s là phân v m c 1 − α/2 c a phân như sau: [ x - t α ], trong ñó t α ; x + tα , n −1 , n −1 , n −1 n n 2 2 2 ph i Student v i b c t do n −1. 4
  5. 2. T ch c ñ T n s xu t hi n c a s li u trong các kho ng cách ñ u nhau cho phép phác ho bi u ñ t n s , còn g i là t ch c ñ . ð v t ch c ñ c n ph i ti n hành phân t / nhóm s l i u. 2.1. T o mi n phân t ð ti n hành phân t s li u (t o Bin), c n th c hi n các bư c sau: - Dùng các hàm Min, Max ñ xác ñ nh giá tr nh nh t và giá tr l n nh t. - ð nh ra giá tr c n dư i và giá tr c n trên c a mi n phân t . - Ghi giá tr c n dư i vào ô ñ u c a mi n phân t và bôi ñen toàn mi n này. - Ch n Edit > Fill > Series ñ khai báo các m c: + Trong m c Series in ch n Columns ( d li u theo c t) + Trong m c Type ch n Linear ( d li u tăng theo c p s c ng) + Trong Step value: nh p giá tr bư c tăng + Trong Stop value: nh p giá tr c n trên + OK. Ví d 2: D a trên 30 s li u v chi u dài cá ta t o mi n phân t (Bin) như trên hình I.3 v i mi n phân t t ô D2 t i ô D12 (k c nhãn), giá tr c n dư i là 10, c n trên là 55, giá tr bư c tăng 5. Hình I.3. T o mi n Bin cho các s li u v chi u dài cá 2.2. V t ch c ñ a. Các bư c th c hi n Ch n Tools> Data Analysis> Histogram ñ khai báo các m c: - Input range: mi n d li u. - Input Bin: mi n phân t . - Labels : nhãn hàng ñ u n u có. - Output range: Mi n k t qu . 5
  6. - Pareto: t n s s p x p gi m d n. - Cumulative Percentage: T n su t c ng d n %. - Chart output: Bi u ñ . - OK. Trong ví d 2 ch n Tools> Data Analysis> Histogram và khai báo như trong hình I.4. Hình I.4. Các khai báo ñ v t ch c ñ b. K t qu v t ch c ñ Hình I.5. T ch c ñ c. Phân tích k t qu - T n s s li u rơi vào t ng kho ng ñư c ghi c n trên c a kho ng. (Ch ng h n, có 2 s li u thu c vào kho ng (10,15], vì v y s 2 ñư c ghi tương ng v i s 15 là c n trên). - Nhìn vào hình I.5. ta có th th y trong kho ng nào s li u xu t hi n nhi u nh t. Ngoài ra, hình d ng c a t ch c ñ còn cho bi t: dãy s li u kh o sát ñư c v chi u dài c a cá có th coi là tuân theo lu t chu n. 6
  7. 3. Tính h s tương quan và tìm phương trình h i qui 3.1. Tính h s tương quan Excel cho phép tính h s tương quan ñơn gi a các bi n s p x p thành m t b ng g m n hàng, n c t (m i c t là 1 bi n). a. Các bư c th c hi n Ch n Tools>Data Analysis>Correlation và khai báo các m c: - Input range: mi n d li u k c nhãn. - Grouped by: Column (s li u theo c t). - Labels in first row : ðánh d u √ vào ô này n u có nhãn hàng ñ u. - Output range: mi n ra. - OK. Ví d 3: ð nghiên c u m i tương quan gi a các ñ c tính dài bông, s h t, s bông v i năng su t lúa, c n th c hi n các bư c sau: - Ch n Tools>Data analysis>Correlation. - Khai báo các m c (xem hình I.6). Hình I.6. Các bư c khai báo khi tính h s tương quan - K t qu thu ñư c trên b ng I.1. B ng I.1. K t qu tính h s tương quan Năng su t Dài bông P1000 S bông Dài bông 1 P1000 0.233314 1 S bông -0.22056 0.340772 1 Năng su t 0.200805 0.66632 0.661379 1 7
  8. b. Phân tích k t qu - H s tương quan c a hàng và c t ghi ô giao gi a hàng và c t. - H s tương quan âm ( < 0) th hi n m i tương quan ngh ch bi n (ch ng h n tương quan gi a “dài bông” và “s bông” là ngh ch bi n). - Các h s tương quan có giá tr tuy t ñ i x p x 0.75 tr lên th hi n m i tương quan tuy n tính m nh gi a hai bi n (tương quan gi a “năng su t” và “P1000” có th t m coi là tương quan tuy n tính m nh). 3.2. Tìm phương trình h i quy Excel cho phép tìm phương trình h i quy tuy n tính ñơn y = a+ bx và h i quy tuy n tính b i y = a0 + a1x1 + a2x2 + . . . + anxn. Các bi n ñ c l p ch a trong n c t, bi n ph thu c y ñ trong m t c t, các giá tr tương ng gi a bi n ñ c l p và bi n ph thu c ñư c x p trên cùng m t hàng. a. Các bư c th c hi n Ch n Tools>Data Analysis>Regression và khai báo các m c: - Input y range: mi n d li u bi n y. - Input x range: mi n d li u các bi n x. - Label: ðánh d u √ vào ô này n u có nhãn hàng ñ u. - Confidence level : 95% ( ñ tin c y 95%). - Constant in zero: ðánh d u √ n u h s t do a0 = 0 . - Output range: mi n xu t k t qu . - Residuals : ðánh d u √ vào ô này ñ hi n ph n dư hay sai l ch gi a y th c nghi m và y theo h i quy. - Standardized residuals: ðánh d u √ ñ hi n ph n dư ñã chu n hoá. - Residuals plot: ðánh d u √ ñ hi n ñ th ph n dư. - Line fit plots: ðánh d u √ ñ hi n ñ th các ñư ng d báo. - Normal probability plot: ðánh d u √ ñ hi n ñ th ph n dư ñã chu n hoá. - OK. Ví d 4: Tìm phương trình h i qui y= a0 + a1x1 + a2yx2 + a3x3 c a năng su t lúa y ph thu c tuy n tính vào ñ dài bông (x1), tr ng lư ng 1000 h t (x2) và s bông / m t cây (x3) v i các s li u cho trong hình I.7. Ch n Tools>Data Analysis>Regression và khai báo các m c như trên hình I.7. ñ thu ñư c k t qu như trên hình I.8. b. Phân tích k t qu - N u h s tương quan b i x p x 0.75 ho c l n hơn thì mô hình h i quy tuy n tính là thích h p (ngư c l i nên tìm mô hình khác). Trong ví d 4 h s tương quan b i là 0.8375 nên mô hình tuy n tính ñư c coi là thích h p. - H s tương quan R square trong ví d 4 là 0.7014 cho bi t 70.14% s bi n ñ ng c a y là do các y u t x1, x2, x3 gây nên. H s Adjusted R square là 62.00% không sát g n v i R square ch ng t không ph i t t c các bi n ñưa vào là th c s c n thi t. 8
  9. - F th c nghi m là 8.6142 ng v i xác su t 0.00316 nh hơn m c xác su t ý nghĩa 0.05 nên phương trình h i quy tuy n tính ñư c ch p nh n. - Nhìn vào các h s c a các bi n ta vi t ñư c ñư ng h i quy d báo. Trong ví d 4 phương trình h i quy là: y = − 3.61899 + 0.085345x1 + 0.081163x2 + 0.02083x3 . Tuy nhiên căn c vào các xác su t cho c t P-value thì h s c a x1 là không ñáng tin c y, vì xác su t tương ng > 0.05 (m c ý nghĩa ñã ch n). Trong trư ng h p này, c n ti n hành l c b t bi n x1 ñ ñư c ñư ng h i quy v i các h s ñ u có ý nghĩa. Hình I.7. Khai báo ñ tìm phương trình h i quy Hình I.8. K t qu tìm phương trình h i quy 9
  10. 4. Phân tích phương sai Phân tích phương sai là công c ch y u ñ phân tích các s li u khi theo dõi nh hư ng c a các nhân t (factor) trong thí nghi m và nh hư ng tương tác c a chúng lên m t (hay nhi u) ch s ñ u ra. ð thu th p s li u, thí nghi m c n ñư c thi t k phù h p v i m c ñích nghiên c u và ñi u ki n c th nơi ti n hành thí nghi m. ð phân tích m t nhân t , thí nghi m thư ng ñư c thi t k theo ki u hoàn toàn ng u nhiên, ki u kh i hoàn toàn ng u nhiên, hay ô vuông La tinh. ð phân tích hai nhân t , thí nghi m ñư c b trí theo ki u tr c giao, ki u chia ô l n, ô v a, ô nh , ho c k t h p v a chia băng v a chia ô. T ba nhân t tr lên thì c n b trí thí nghi m sao cho m i nhân t có hai m c hay m i nhân t có ba m c. 4.1. Phân tích phương sai m t nhân t Phân tích phương sai m t nhân t ñư c s d ng ñ phân tích s li u khi theo dõi nh hư ng c a các m c c a nhân t t i k t qu , như nh hư ng c a các công th c cho ăn ñ n năng su t th t l n, nh hư ng c a các công th c phun thu c sâu ñ n t l sâu b nh... ð phân tích phương sai m t nhân t c n thi t k thí nghi m ki u hoàn toàn ng u nhiên, m i m c l p l i m t s l n, s l n l p c a các m c c a nhân t không c n ph i b ng nhau. Thi t k thí nghi m m t nhân t hoàn toàn ng u nhiên (CRD). G i k m c c a nhân t hay k công th c c n ti n hành là T1, T2 . . . , Tk . Ch ng h n trong thí nghi m xem xét nh hư ng c a 11 lo i thu c phòng sâu b nh t i năng su t c a m t gi ng lúa, nhân t ñây ch g m m t y u t có 11 m c là 11 lo i thu c nên k = 11. M i lo i thu c ñư c th nghi m trên m t s ô thí nghi m (hay ñơn v thí nghi m), m i ô ñư c coi là m t l n l p. N u thí nghi m 5 gi ng lúa và 11 lo i thu c trên và ch xét tác ñ ng chung c a t h p gi ng và thu c (Gi × Pj) t i năng su t lúa thì có thí nghi m m t nhân t v i k = 5× 11 = 55 công th c thí nghi m. S ô thí nghi m (hay s l n l p) cho m i công th c có th ch n tuỳ ý, không nh t thi t ph i b ng nhau. Phân tích phương sai m t nhân t ñư c ti n hành v i các d li u ñư c s p thành nhi u nhóm, m i nhóm là các l n l p c a m t m c c a nhân t , nh m tách bi t các phương sai theo hai ngu n bi n ñ ng nhân t và sai s . V i i = 1, 2, …, k, m i công th c Ti ñư c th c hi n trên ni ô thí nghi m, các k t qu th c nghi m xij ñư c coi như m t m u th c nghi m ñ i v i bi n ng u nhiên Xi. D a vào k t qu th c nghi m c n ñưa ra suy ñoán v vi c các trung bình mi c a các bi n Xi là như nhau (t c là các công th c không nh hư ng gì ñáng k t i ch s c n kh o sát) hay là khác nhau. Có nhi u ki u thi t k thí nghi m ñ gi i quy t bài toán này. Gi s nhân t có a m c, m c i ñư c l p l i ni l n, như v y t ng s có n = ∑ ni quan sát, hay còn nói là có n ô thí nghi m. N u b trí n ô thí nghi m hoàn toàn ng u nhiên ta có thi t k thí nghi m hoàn toàn ng u nhiên (completely randomized design). Khi ti n hành thí nghi m ki u này ph i dùng n phi u ghi t 1 ñ n n, rút thăm ng u nhiên n1 phi u ñ có các ô thí nghi m ñ i v i công th c 1, rút ti p n2 phi u ñ có các ô thí nghi m ñ i v i công th c 2, ..., nk ô cu i cùng là c a công th c k. Vi c rút thăm ng u nhiên ñư c th c hi n trên toàn b các ô thí nghi m. Vi c tính toán và k t lu n d a trên mô hình: xij = µ + αi + eij (i = 1, …, k và j = 1, ..., ni), v i xij là k t qu c a l n l p th j c a m c i, µ là trung bình chung, αi là nh hư ng c a m c i c a nhân t , còn eij là sai s ng u nhiên. xij có trung bình mi = µ +αi. Các sai s 10
  11. eij ñư c gi thi t là ñ c l p và tuân theo phân ph i chu n v i kỳ v ng 0 và phương sai σ2. k ∑α Các αi ñư c coi là tho mãn ñi u ki n = 0. i i =1 a. Các bư c th c hi n S li u ñư c ñi n theo c t ho c theo hàng (n u vào theo hàng thì m i hàng ng v i m t m c c a nhân t ), ô ñ u tiên ghi tên m c, các ô ti p theo ghi s li u. Ch n Tools> Data Analysis > Anova: Single Factor và khai báo: - Input range: Khai báo mi n d li u vào (m t ch nh t bao trùm toàn b các ô ch a tên m c và toàn b các s li u). - Grouped by: Column (s li u theo c t) ho c row (s li u theo hàng). - Label in First column : nhãn hàng ñ u. - Alpha: 0.05 (m c ý nghĩa α). - Output range: mi n ra. b. Phân tích k t qu - K t qu in ra g m các th ng kê cơ b n cho t ng m c (trung bình, ñ l ch chu n...) và b ng phân tích phương sai. - N u giá tr xác su t P-value < alpha (ho c F th c nghi m > F lý thuy t) thì các công th c có tác ñ ng khác nhau t i k t qu , ngư c l i các công th c không có khác bi t ñáng k . - N u k t lu n các công th c có tác ñ ng khác nhau t i k t qu thì ph i ti n hành bư c ti p theo là so sánh các công th c ñ rút ra công th c nào t t nh t. Ví d 5: Thí nghi m nh hư ng c a các lo i thu c ñ n năng su t lúa (11 lo i thu c là T1 ñ n T11, 4 c t s li u là năng su t thu ñư c), s li u thu ñư c cho trong b ng I.2., các l nh th c hi n trong Excel ñư c minh ho trong hình I.9., còn k t qu cho trong b ng I.3. B ng I.2. nh hư ng c a các lo i thu c ñ n năng su t (ns) lúa Lo i thu c NS ô 1 NS ô 2 NS ô 3 NS ô 4 T1 3.187 4.61 3.562 3.217 T2 3.39 2.875 2.775 T3 2.797 3.001 2.505 3.49 T4 2.832 3.103 3.448 2.255 T5 2.233 2.743 2.727 T6 2.952 2.272 2.47 T7 2.858 2.895 2.458 1.723 T8 2.308 2.335 1.957 T9 2.013 1.788 2.248 2.115 T10 3.202 3.06 2.24 2.69 T11 1.192 1.652 1.075 1.03 OK 11
  12. Hình I.9. Th c hi n phân tích phương sai m t nhân t trong Excel. B ng I.3. K t qu phân tích phương sai Groups Count Sum Average Variance T1 4 14.576 3.644 0.443686 T2 3 9.04 3.013333 0.108908 T3 4 11.793 2.94825 0.171874 T4 4 11.638 2.9095 0.253934 T5 3 7.703 2.567667 0.084065 T6 3 7.694 2.564667 0.122321 T7 4 9.934 2.4835 0.296198 T8 3 6.6 2.2 0.044469 T9 4 8.164 2.041 0.037706 T 10 4 11.192 2.798 0.184963 T 11 4 4.949 1.23725 0.081114 ANOVA Source of SS df MS F P-value F crit Variation Between 15.1039 10 1.51039 8.54171 2.66E-06 2.1768 Groups Within 5.1279 29 0.17682 Groups Total 20.2319 39 12
  13. T b ng I.3 ta k t lu n các công th c có tác ñ ng khác nhau t i năng su t lúa. Gi i thích: Như ñã nói trên, phân tích phương sai m t nhân t tách bi t các phương sai theo hai ngu n bi n ñ ng nhân t và sai s . Theo b ng I.3, ta có: - T ng sai l ch gi a các nhóm do các m c c a nhân t t o nên là SST (sum of 2 2  ni   k ni  ( ) k k 2 squares treatment) = ∑ n i x i − x = ∑  ∑ x ij  / n i −  ∑∑ x ij  / n = 15.1039. Do ñó i =1  j=1   i=1 j=1  i =1 t ng sai l ch trung bình gi a các nhóm là MSST = SST/ (k−1) = 15.1039/ 10 = 1.5139. - T ng sai l ch trong t ng nhóm do các sai s eij gây nên là SSE (sum of squares 2  k ni 2  k  ni  ni error) = ∑∑ ( x ij − x ) =  ∑∑ x ij  − ∑  ∑ x ij  / n i = 5.1279. Do ñó t ng sai l ch trung k 2  i=1 j=1  i =1  j=1  i =1 j=1 bình trong t ng nhóm là MSSE = SSE/ (n−k) = 5.11279/ 29 = 0.12682. V i các gi thi t ñã nêu, có th ch ng minh ñư c ñ i lư ng th ng kê F = MSSB/MSSE (ñ i v i m u lý thuy t) tuân theo phân ph i Fisher v i b c t do là (k−1, n−k). Do ñó, n u F th c nghi m = 8.54171 > F lý thuy t = f (0.05; 10, 29) = 2.1768 thì gi thuy t “các công th c s d ng thu c không nh hư ng t i năng su t (trung bình) c a lúa” m c ý nghĩa α = 0.05. ði u này có nghĩa là các công th c s d ng thu c có tác b bác b ñ ng khác nhau t i năng su t lúa. Ngoài ra, t b ng I.3 có th nh n xét r ng công th c T1 cho năng su t cao nh t.Tuy nhiên, ñ k t lu n chính xác hơn v ñi u này c n so sánh trung bình gi a các nhóm. B ng I.4. So sánh trung bình gi a các nhóm T1-Ti Groups Count Sum Average T1 4 14.576 3.6440 T2 3 9.040 3.0133 0.6307 T3 4 11.793 2.9483 0.6958 T4 4 11.638 2.9095 0.7345 T10 4 11.192 2.7980 0.8460 T5 3 7.703 2.5677 1.0763 T6 3 7.694 2.5647 1.0793 T7 4 9.934 2.4835 1.1605 T8 3 6.600 2.2000 1.4440 T9 4 8.164 2.0410 1.6030 T11 4 4.949 1.2373 2.4068 Ph n m m Excel không cho phép so sánh các trung bình c a các nhóm ng v i các m c c a nhân t (các công th c). Tuy nhiên, n u c n so sánh trung bình mi (v i ni l n l p) v i trung bình mj (nj l n l p) ngư i dùng có th t tính thêm LSD (Least Significance Difference) theo công th c LSD = tα,df × SQRT(s2(1/ni + 1/nj)), trong ñó s2 là phương sai chung ñư c ư c lư ng b i trung bình sai s bình phương trong n i b nhóm (MS within 13
  14. groups), α = 1-p, và tα , df là giá tr t c a b ng Student ng v i m c ý nghĩa α và df b c t do. tα , df có th tìm ñư c b ng cách tra b ng s hay b ng hàm TINV trong Excel. Trong ví d 5, ñ so sánh nh hư ng c a thu c T1, T2 ñ n năng su t lúa, trư c h t c n tính tr tuy t ñ i | m1- m2| c a hi u các năng su t trung bình m1, m2 khi s d ng 2 lo i thu c trên (chính b ng 0.6307, xem b ng I.4). ng v i t = t(0.05 , 29) = 2.045 (tra t b ng Student)có s2= 0.17682. LSD ñư c tính cho các trư ng h p ri, rj b ng 3 ho c 4 như sau: LSD= 2.045 × SQRT( 0.17682×(1/3+1/4) = 0.656739049; LSD= 2.045 × SQRT( 0.17682×(1/4+1/4) = 0.608022212; LSD= 2.045 × SQRT( 0.17682×(1/3+1/3) = 0.702083575. Trong trư ng h p tính nh hư ng c a thu c T1, T2 t i năng su t trung bình c a lúa ta có: | m1- m2| = 0.6307 <LSD = 0.656739049 nên công th c T1, T2 không khác nhau rõ r t. Các k t qu tính toán tương t ñư c ghi trong b ng I.3. cho th y công th c T1 và T2 là không khác nhau rõ r t, công th c T1 khác các công th c t T3 ñ n T11. Tương t có th so sánh công th c T2 v i các công th c t T3 ñ n T11. Công th c T1 cho năng su t cao nh t là t t nh t, công th c T11 cho năng su t nh nh t là kém nh t. 4.2. Phân tích phương sai hai nhân t không tương tác Khi phân tích phương sai hai nhân t A và B có th x y ra các hai trư ng h p: trư ng h p A và B không tương tác (bi n ñ ng gây nên b i tác ñ ng ñ ng th i c a A và B g n sát 0) và trư ng h p A và B tương tác (n u trái l i). Phân tích phương sai m t nhân t b trí ki u kh i hoàn toàn ng u nhiên ñư c coi là trư ng h p riêng c a phân tích phương sai hai nhân t không tương tác (nhân t kh i là nhân t th hai không tương tác v i nhân t th nh t). Thi t k thí nghi m theo kh i ng u nhiên ñ y ñ RCBD. Gi s có k công th c, m i công th c l p l i r l n. T t c có n = k×r ô thí nghi m. ð tránh các tác ñ ng c a m t s y u t ngo i c nh lên k t qu ñ u ra c a m t s công th c nào ñó, chúng ta c n b trí các ô thí nghi m m t cách h p lý. Trong ñi u ki n không có ñ n ô thí nghi m ñ ng ñ u, c n thi t k thí nghi m b ng cách chia thí nghi m thành r kh i v i k ô trong m i kh i tương ñ i ñ ng ñ u v m i m t sao cho các tác ñ ng ph không nh hư ng t i thí nghi m. Sau ñó, xét kh i th nh t và làm k phi u ñ b t thăm xem k công th c x p vào k ô nào. Ti p t c b t thăm cho kh i th hai, th ba, . . . cho t i kh i th r. Vi c chia kh i thí nghi m nói chung ph thu c vào ñ a ñi m thí nghi m. Ch ng h n, c n chia kh i th ng góc v i m t hư ng bi n ñ ng có nh hư ng ñ n k t qu thí nghi m như hư ng gió, hư ng ch y c a nư c ng m, hư ng n ng, hư ng d c, hư ng thay ñ i c a ñ phì c a ñ t ... sao cho m i công th c có m t m t l n m t m c c a bi n ñ ng. Vi c chia kh i thí nghi m cũng có th ph thu c vào th i gian ti n hành thí nghi m v i các tác ñ ng c a th i ti t. N u m i ngày ch làm ñư c k thí nghi m và ta ch có r ngày ñ làm t t c các thí nghi m, thì c n ph i phân chia vi c th c hi n các thí nghi m ra r ngày, như v y ñây ngày là kh i. M t cái l i n a là trong thi t k thí nghi m RCBD có th ch n kh i khác nhau v không gian ho c khác nhau v th i gian (nhưng không ñư c khác nhau quá xa ñ n m c có s thay ñ i ñi u ki n thí nghi m). Do ñó, k t lu n rút ra có tính khái quát cao hơn khi so v i k t lu n ñ t ñư c trong thi t k thí nghi m hoàn toàn ng u nhiên (t p trung toàn b các thí nghi m vào m t nơi hay cùng m t th i gian). 14
  15. Vi c tính toán và k t lu n d a trên mô hình: xij = µ + αi + βj + eij (i = 1, …, k và j = 1, ..., r), v i xij là k t qu c a m c i kh i j, µ là trung bình chung, αi là nh hư ng c a m c i c a nhân t , βj là nh hư ng c a kh i j, còn eij là sai s ng u nhiên. Các sai s eij ñư c gi thi t là ñ c l p và tuân theo phân ph i chu n v i kỳ v ng 0 và phương sai σ2. Các k r ∑ αi = ∑β tham s αi và βj ñư c coi là tho mãn ñi u ki n = 0. j i =1 j=1 a. Các bư c th c hi n Khi phân tích phương sai hai nhân t không tương tác, s li u c n ñư c s p x p theo cách sau: hàng là các m c c a nhân t th nh t, c t là các m c c a nhân t th hai (trong trư ng h p c n phân tích phương sai m t nhân t b trí ki u kh i ng u nhiên thì hàng là các m c c a nhân t , c t là các kh i ng u nhiên). Ch n Tools >Data Analysis >Anova: Two Factor Without Replication sau ñó khai báo ti p các thông tin trong hình I.10 và kích OK. Hình I.10. H p tho i khai báo ñ phân tích phương sai không tương tác Ví d 6: B trí thí nghi m phân tích nhân t 1 (có b n m c) theo kh i hoàn toàn ng u nhiên (nhân t 2 có b n kh i), ta thu ñư c các s li u như trong b ng I.5. B ng I.5. S li u phân tích m t nhân t theo kh i ng u nhiên Kh i 1 Kh i 2 Kh i 3 Kh i 4 M c1 47 52 62 51 M c2 50 54 67 57 M c3 57 53 69 57 M c4 54 65 74 59 K t qu thu ñư c khi phân tích phương sai cho b ng I.6. 15
  16. b. Phân tích k t qu - Các m c c a nhân t 1 có nh hư ng khác nhau ñ n k t qu (F th c nghi m > F lý thuy t). - Các m c c a nhân t 2 có nh hư ng khác nhau ñ n k t qu (F th c nghi m > F lý thuy t). Chú ý: F lý thuy t có th tìm b ng hàm FINV(0.05, 3, 9) =3.86254, và giá tr t cũng có th tìm ñư c b ng hàmTINV(0.05, 9) =2.262. B ng I.6. K t qu phân tích phương sai hai nhân t không tương tác Anova: Two-Factor Without Replication Phân tích nhân t 1(hàng) SUMMARY Count Sum Average Variance M c1 4 212 53 40.6667 M c2 4 228 57 52.6667 M c3 4 236 59 48 M c4 4 252 63 74 Phân tích nhân t 2 (c t) Kh i 1 4 208 52 19.3333 Kh i 2 4 224 56 36.6667 Kh i 3 4 272 68 24.6667 Kh i 4 4 224 56 12 ANOVA Source of SS df MS F P-value F crit Variation Rows 208 3 69.3333 8.91429 0.00465 3.86254 Columns 576 3 192 24.6857 0.00011 3.86254 Error 70 9 7.77778 Total 854 15 Gi i thích: Phân tích phương sai hai nhân t không tương tác tách bi t các phương sai theo ba ngu n bi n ñ ng nhân t A (kh i), nhân t B và sai s . Theo b ng I.6, ta có: - T ng sai l ch toàn ph n là SSTO (total sum of squares) = ∑∑ ( x − x ) = ∑∑ x ij − nx 2 = 854. k r r k 2 2 ij i =1 j=1 j=1 i=1 - T ng sai l ch do nhân t A (kh i) là SSB (sum of squares due to block) 2 ( ) 1 rk  r 2 ∑  ∑ x ij  − nx 2 = 576. Do ñó t ng sai l ch trung bình gi a các nhóm r ∑ x. j − x = k j=1  i=1  j=1 c a nhân t A là MSSB = SSB/ (r−1) = 576/ 3 = 132. 16
  17. - T ng sai l ch do nhân t B là SST (sum of squares due to treatment) 2 1kr  ( ) k 2 k ∑ x i. − x ∑  ∑ x ij  − nx 2 = 208. Do ñó t ng sai l ch trung bình gi a các nhóm = r i=1  j=1  i =1 c a nhân t B là MSST = SST/ (k−1) = 208 / 3 = 69.3333. - T ng sai l ch do sai s ng u nhiên là SSE (sum of squares due to errors) = SSTO − SSB − SST = 854 − 576 − 208 = 70. Do ñó t ng sai l ch trung bình c a sai s ng u nhiên là MSSE = SSE/ (n−k−r+1) = 70 / 9 = 7.7778. V i gi thi t ñã nêu, có th ch ng minh ñư c r ng ñ i lư ng th ng kê FR = MSST/MSSE (ñ i v i m u lý thuy t) tuân theo phân ph i Fisher v i b c t do là (k−1, n−k−l+1). Do ñó, n u FR th c nghi m = 8.91429 > FR lý thuy t = f (0.05; 3, 9) = 3.86254 thì gi thuy t “các công th c không nh hư ng t i s trung bình chung” không ñư c ch p nh n theo quy t c ki m ñ nh có m c ý nghĩa α = 0.05. ði u này có nghĩa là các công th c khác nhau có tác ñ ng khác nhau t i s trung bình chung. Ngoài ra, c n chú ý r ng do FB th c nghi m = 24.6857 > FB lý thuy t = f (0.05; 3, 9) = 3.86254 nên gi thuy t “các công th c không nh hư ng t i s trung bình chung” cũng b bác b theo quy t c ki m ñ nh có m c ý nghĩa α = 0.05. ði u này có nghĩa là các kh i khác nhau có tác ñ ng khác nhau t i s trung bình chung. Chú ý: Vi c thi t k thí nghi m kh i ng u nhiên ñ y ñ là r t h p lý khi ch g p m t y u t h n ch , t c là ch có m t ngu n bi n ñ ng duy nh t nh hư ng t i khu v c thí nghi m. Trong trư ng h p có hai y u t h n ch (hai ngu n bi n ñ ng) nh hư ng t i thí nghi m, thí nghi m có th ñư c thi t k theo ki u ô vuông La tinh (Latin squares). Các s li u thu th p ñư c thu th p theo hàng (m i hàng là m t m c c a nhân t A, ch ng h n như hư ng gió trong thí nghi m kh o sát năng su t các lo i lúa ch u h n) và theo c t (m i c t là m t m c c a nhân t B, ch ng h n như ñ cao c a ñ a ñi m canh tác). Trong thi t k thí nghi m theo ô vuông La tinh, s m c c a nhân t A b t bu c b ng s m c c a nhân t B và b ng k. Ngoài ra, s công th c c n ti n hành thí nghi m (s m c c a nhân t C) cũng b ng k. Thí nghi m ñư c thi t k ng u nhiên sao cho t i m i t h p (m t m c c a nhân t A, m t m c c a nhân t B) có duy nh t m t công th c thí nghi m (m t m c c a nhân t C) ñư c ti n hành. Như v y thay vì s ô thí nghi m là k3, trong thi t k thí nghi m ki u ô vuông La tinh chúng ta ch c n có k2 ô thí nghi m. Các công th c tính SSTO (t ng sai l ch toàn ph n), SSA , SSB, SSC và SSE (các t ng sai l ch do các nhân t A, B, C và sai s ng u nhiên) cũng ñư c tính tương t như các công th c ñã bi t trên ñây. Các t ng sai l ch này có các b c t do tương ng là k2 −1, k −1, k − 1, k −1 và k2 − 3k +2 (v i k ≥ 4). T ñó tính ñư c MSSA, MSSB, MSSC và MSSE. Ti p theo c n thi t l p các giá tr F th c nghi m là: FA = MSSA/MSSE, FB = MSSB/MSSE, FC = MSSC/MSSE ñ rút ra các suy ñoán th ng kê tương ng. S d ng ch c năng phân tích s li u c a Excel cũng có th giúp x lý ñư c các s li u thu ñư c khi thi t k thí nghi m theo ki u ô vuông La tinh m t cách nhanh chóng hơn. Trư c h t c n nh p s li u c a các ô vuông La tinh vào b ng tính Excel (gi s các s li u thu ñư c cúng gi ng như trong ví d 6), sau ñó th c hi n l nh Tools >Data Analysis >Anova: Two Factor Without Replication ñ thu ñư c k t qu tương t như trên b ng I.6. Lúc ñó s có SSA = 208, SSB = 576. Sau ñó c n tính tr c ti p ñ thu ñư c SSC theo công 2 2 1k  1k  x2 th c: ∑  ∑ x ijs  − nx = ∑  ∑ x ijs  − 2 , trong ñó xijs là các s li u thu ñư c khi 2 k s=1  i,j k s=1  i,j k  th c hi n công th c s, v i s = 1, 2, …, k. Gi s các s li u ñư c t ng h p như t i b ng I.7. 17
  18. B ng I.7. S li u phân tích m t nhân t theo ô vuông La tinh Các m c B1 B2 B3 B4 A1 47 (C3) 52 (C4) 62 (C1) 51 (C2) A2 50 (C2) 54 (C3) 67 (C4) 57 (C1) A3 57 (C1) 53 (C2) 69 (C3) 57 (C4) A4 54 (C4) 65 (C1) 74 (C2) 59 (C3) Lúc ñó, SSC = (2412 + 2282 + 2292 + 2302)/4 −(241 + 228 + 229 + 230)2/16 = 27.5. Do ñó, SSE = SSTO − SSA − SSB − SSC = 854 − 208 − 576 − 27.5 = 42.5. T ñó tính ñư c MSSC = 27.5/ 3 = 9.1667 và MSSE = 42.5/ 6 = 7.0834. V y FC th c nghi m = MSSC / MSSE = 7.0834 / 9.1667 = 0.77273. Trong khi ñó FC lý thuy t = f(0.05; 3, 6) = 4.757. Suy ñoán th ng kê có th ñư c ñưa ra là: các công th c không nh hư ng ñáng k t i s trung bình chung c a ch s kh o sát. 4.3. Phân tích phương sai hai nhân t Trong trư ng h p này khi phân tích phương sai, ngoài tác ñ ng c a t ng nhân t A và nhân t B lên k t qu c a thí nghi m, ta ph i tính ñ n s tác ñ ng ñ ng th i còn g i là tác ñ ng tương tác c a c hai nhân t này. Thi t k thí nghi m hai nhân t . M t s ki u thi t k thí nghi m ñư c áp d ng ñ phân tích phương sai hai nhân t tương tác là: thi t k thí nghi m tr c giao (hai nhân t chéo nhau , OED), thi t k thí nghi m phân c p (hai nhân t l ng nhau, HED), thi t k thí nghi m chia ô (SPED) và thi t k thí nghi m chia băng. Hình I.11 minh ho các cách thi t k thí nghi m. Tuỳ theo m c ñích và ñi u ki n thí nghi m trong các lĩnh v c chuyên môn, thí nghi m ñư c thi t k theo cách th c thích h p và s li u th c nghi m thu ñư c cũng ñư c x lý m t cách phù h p nh m rút ra các suy ñoán th ng kê có ý nghĩa. 4.3.1 Thi t k thí nghi m tr c giao Trư ng h p ñơn gi n nh t c a mô hình chéo nhau là y u t A có 2 m c A1 và A2, y u t B có 2 m c B1 và B2. Các t h p có th c a các m c y u t là: Y ut B Y ut A B1 B2 A1 A1B1 A1B2 A2 A2B1 A2B2 4.3.2 Thi t k thí nghi m phân c p Ki u thí nghi m hai nhân t phân c p (Hierachical) hay chia (Nested) thư ng ñư c dùng trong các nghiên c u v di truy n. Trong ñó m t nhân t là c p trên, m t nhân t là c p dư i, thí nghi m l p l i r l n. ð c th xét thí d A là bò ñ c gi ng, t t c có 4 con A1, A2, A3, A4. M i con ñ c cho ph i v i 3 con cái g i t t là B1, B2, B3. M i con bò cái sinh 4 con. Ta có sơ ñ sau: 18
  19. A 1 2 3 4 B 1 2 3 1 2 3 1 2 3 1 2 3 x111 x121 x131 x211 x221 x131 x311 x321 x331 x411 x421 x431 x112 x122 x132 x212 x222 x232 x312 x322 x332 x412 x422 x432 x113 x123 x133 x213 x223 x233 x313 x323 x333 x413 x423 x433 x114 x124 x134 x214 x224 x234 x314 x324 x334 x414 x424 x434 4.3.2 Thi t k thí nghi m phân c p theo kh i Thư ng b trí thí nghi m theo kh i, m i kh i chia thành a ô l n ñ b t thăm cho a m c c a nhân t A. Vi c b t thăm ñư c th c hi n riêng r cho t ng kh i. M i ô l n chia thành b ô nh ñ b t thăm cho b m c c a nhân t B. Vi c b t thăm th c hi n riêng r cho t ng ô l n. Thí d y u t A có 4 m c (A1, A2, A3và A4), y u t B có 2 m c (B1 và B2). Ba m c c a y u t A ñư c b trí trên ô l n trong 3 kh i. M i ô l n chia nh thành 2 ô nh ñ b trí ng u nhiên các m c c a y u t B. Sơ b trí thí nghi m có th ñư c trình bày như sau: Kh i 1 Kh i 2 Kh i 3 A4 A1 A2 A3 A2 A1 A4 A3 A1 A2 A4 A3 B2 B2 B1 B2 B1 B2 B1 B1 B2 B1 B2 B1 B1 B1 B2 B1 B2 B1 B2 B2 B1 B2 B1 B2 Thi t k thí nghi m tr c giao. Chúng ta ñi sâu vào thi t k thí nghi m tr c giao. Gi s nhân t A có k m c là A1, A2, …, Ak và nhân t B có r m c là B1, B2, … Br. S công th c là k×r, m i công th c ñư c l p l i s l n. Như v y chúng ta có t t c k×r×s ô thí nghi m. Có th thi t k thí nghi m tr c giao theo ki u ng u nhiên hoàn toàn (CRD) ho c theo ki u kh i ng u nhiên ñ y ñ (RCBD). Trong trư ng h p th nh t ta c n b t thăm các ô thí nghi m ñ phân vào m i ô m t công th c: trư c h t b t thăm ng u nhiên s ô ñ phân công cho công th c th nh t, ti p theo b t thăm s ô ñ phân cho công th c th c 2, …, làm như v y cho t i công th c th k×r. Trong trư ng h p th hai, ta c n b trí ñ s kh i, m i kh i ph i có ñ k×r công th c ñư c phân vào các ô m t cách ng u nhiên. Vi c tính toán và k t lu n d a trên mô hình: xijq = µ + αi + βj + (αβ)ij + eijq (i = 1, …, k, j = 1, ..., r và q = 1, 2, …, s), v i xijq là k t qu c a các m c i c a nhân t A, m c j c a nhân t B và ô thí nghi m th q, µ là trung bình chung, αi là nh hư ng c a m c i c a nhân t A, βj là nh hư ng m c j c a nhân t B, (αβ)ij là nh hư ng c a s tương tác c a m c i c a A và m c j c a B, còn eijq là sai s ng u nhiên. Các sai s eijq ñư c gi thi t là ñ c l p và tuân theo phân ph i chu n v i kỳ v ng 0 và phương sai σ2. Các tham s αi và βj k k r r ∑ α = ∑ β = ∑ ( αβ ) ∑ ( αβ ) ñư c coi là tho mãn ñi u ki n = =0. i j ij ij i =1 j=1 i=1 j=1 19
  20. Sau khi ti n hành thí ngh m, s li u thu ñư c ñư c s p x p như sau: - Nhân t A ñánh theo hàng v i các m c khác nhau, nhân t B ñánh theo c t v i các m c khác nhau. - M i m c c a nhân t A ñư c dành s hàng (cho s l n l p) còn m i m c c a nhân t B ñư c dành ñúng 1 c t. - Tên m i m c c a nhân t A ch vi t m t l n trong s ô c t ñ u, còn tên các m c c a nhân t B thì ghi ñ u m i c t trên hàng ñ u, k t c t th 2. Các ô t hàng 2 c t 2 tr ñi ghi k t qu c a các l n l p c a các t t h p m c. - B ng phân tích phương sai có 5 hàng: Hàng cho nhân t A, hàng dành cho nhân t B, hàng cho tương tác A× B, hàng cho sai s và hàng Total. Giá tr F lý thuy t ñư c tính c t cu i, ta có th ki m tra qua hàm FINV (trong Excel). C t P – value là xác su t tương ng v i giá tr F th c nghi m, n u giá tr này nh hơn alpha thì ta k t lu n nhân t (ho c tương tác) tương ng có nh hư ng ñ n k t qu thí nghi m. Mu n so sánh các trung bình ta làm như phân tích m t nhân t sau khi tìm giá tr t b ng hàm TINV v i s b c t do c a sai s và căn c vào s l n l p c a các trung bình mà ta mu n so sánh. a. Các bư c th c hi n Ch n Tools >Data Analysis >Anova: Two Factor With Replication, sau ñó khai báo các thông tin như trong hình I.12 và kích OK. Ví d 7: Nghiên c u nh hư ng c a vi c bón phân khoáng (nhân t A) theo b n công th c và m t ñ tr ng (nhân t B) g m ba m c t i s n lư ng bông ta có b ng s li u (b ng I.8). B ng I.8. Năng su t bông (t /ha) M c1 M c2 M c3 C.th c 1 14 15 19 15 17 19 16 19 18 21 18 17 C.th c 2 20 18 21 19 19 20 23 18 21 19 20 23 C.th c 3 21 21 21 19 22 22 22 21 18 20 23 21 C.th c 4 20 21 24 23 22 23 21 19 21 19 20 25 20
Theo dõi chúng tôi
Đồng bộ tài khoản