Какие выборки считаются связанными Опции

[Stefa]

Всем добрый день! Столкнулась с проблемой и просто зациклилась на ней. Пример из книги В.Ю. Урбаха : две делянки пшеницы, одна опыт, вторая контроль, измерялась урожайность раз в год.
Год| 1947|1948|1949|1950|1951|1952|1953
Опыт|22.9|20.2|19.5|30.5|35.6|31.9|27.7
Контроль|19.4|16.2|16.9|29.3|31.4|28.5|26.6
Для сравнения урожайности применяется критерий Стьюдента для парных выборок. У меня аналогичная задача, но я не могу доказать, почему эти выборки следует считать парными. На все мои объяснения, что опыт и контроль связаны годом, и что нельзя сравнивать урожайность первой делянки за 1947 год и урожайность второй делянки, например, за 1953г., а необходимо рассматривать именно пары, мне рассказывают про пациентов до и после лечения, и что там да, связанные, а здесь никакой связи нет. Может быть я не права? А если права, то, как объяснить так, чтобы не у кого не возникало никаких сомнений. Помогите, пожалуйста! Заранее большое спасибо.

[Stefa]

Прошу прощения. У меня два города и показатели числа коек на 1000 населения в них за 5 лет:
Год|2006|2007|2008|2009|2010
Город 1|8.5|8.8|7.9|7.5|8.1
Город 2|10.6|12.4|10.9|10.6|7.7

[p2004r]

Прошу прощения. У меня два города и показатели числа коек на 1000 населения в них за 5 лет:
Год|2006|2007|2008|2009|2010
Город 1|8.5|8.8|7.9|7.5|8.1
Город 2|10.6|12.4|10.9|10.6|7.7

1 по моему сравнение с вегетацией растений и зависимостью от года не очень корректно в данном случае

2 модель данных другая, там все в течении года появляется и убирается. в случае койкомест большая их часть переходит от года к году. (может надо каким то образом отображать движение койкомест (прибыло - убыло))?

3 ну поскольку "глобус выдан" попробуем посмотреть на данные

Код

# считываем "как есть"

> read.table("data.txt",sep="|")
       V1     V2     V3     V4     V5     V6
1     Год 2006.0 2007.0 2008.0 2009.0 2010.0
2 Город 1    8.5    8.8    7.9    7.5    8.1
3 Город 2   10.6   12.4   10.9   10.6    7.7


# преобразуем в "широкий" датафрейм

data<-as.data.frame(t(read.table("data.txt",sep="|")[,2:6]))
names(data)<-read.table("data.txt",sep="|")[,1]

> data
    Год Город 1 Город 2
V2 2006     8.5    10.6
V3 2007     8.8    12.4
V4 2008     7.9    10.9
V5 2009     7.5    10.6
V6 2010     8.1     7.7

# преобразуем в "длинный" датафрейм
library(reshape)

> data.long<-melt(data=data,id.vars="Год", measure.vars=c("Город 1","Город 2"))
> data.long
    Год variable value
1  2006  Город 1   8.5
2  2007  Город 1   8.8
3  2008  Город 1   7.9
4  2009  Город 1   7.5
5  2010  Город 1   8.1
6  2006  Город 2  10.6
7  2007  Город 2  12.4
8  2008  Город 2  10.9
9  2009  Город 2  10.6
10 2010  Город 2   7.7

# строим модель смешанных эффектов

library(lme4)

# полную

> model.ful <- lmer(value ~ 1 + (1|Год) + (1|variable), data=data.long)
> model.ful
Linear mixed model fit by REML
Formula: value ~ 1 + (1 | Год) + (1 | variable)
   Data: data.long
   AIC   BIC logLik deviance REMLdev
41.91 43.12 -16.96       36   33.91
Random effects:
Groups   Name        Variance Std.Dev.
Год      (Intercept) 0.3120   0.55857
variable (Intercept) 2.3455   1.53150
Residual             1.2685   1.12628
Number of obs: 10, groups: Год, 5; variable, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept)    9.300      1.167   7.969

# только год

> model.год <- lmer(value ~ 1 + (1|Год) , data=data.long)
> model.год
Linear mixed model fit by REML
Formula: value ~ 1 + (1 | Год)
   Data: data.long
   AIC   BIC logLik deviance REMLdev
43.27 44.17 -18.63    37.79   37.27
Random effects:
Groups   Name        Variance Std.Dev.
Год      (Intercept) 0.0000   0.0000  
Residual             2.8489   1.6879  
Number of obs: 10, groups: Год, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)   9.3000     0.5337   17.42

#только город

> model.город <- lmer(value ~ 1 +  (1|variable), data=data.long)
> model.город
Linear mixed model fit by REML
Formula: value ~ 1 + (1 | variable)
   Data: data.long
   AIC   BIC logLik deviance REMLdev
40.07 40.98 -17.04    36.12   34.07
Random effects:
Groups   Name        Variance Std.Dev.
variable (Intercept) 2.2831   1.5110  
Residual             1.5805   1.2572  
Number of obs: 10, groups: variable, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept)     9.30       1.14   8.158

# отображаем все графически с доверительными интервалами

dotplot(ranef(model.ful, data=data.long, postVar = TRUE))$Год
dotplot(ranef(model.ful, data=data.long, postVar = TRUE))$variable
dotplot(ranef(model.год, data=data.long, postVar = TRUE))
dotplot(ranef(model.город, data=data.long, postVar = TRUE))

# ну и сами данные тоже отображаем графически

print(dotplot(reorder(Год, value) ~ value, data.long, groups = variable, ylab = "Год"))

судя по всему год не имеет достоверного вклада в модель (в качестве фактора дизайна эксперимента), отмечается слабая тенденция.

Книга (на простом английском языке) с описание методик анализа лежит здесь
http://lme4.r-forge.r-project.org/book/

отличия между моделями

Код

> anova(model.ful,model.город)
Data: data.long
Models:
model.город: value ~ 1 + (1 | variable)
model.ful: value ~ 1 + (1 | Год) + (1 | variable)
            Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)
model.город  3 42.116 43.024 -18.058                        
model.ful    4 44.004 45.215 -18.002 0.1121      1     0.7378
> anova(model.ful,model.год)
Data: data.long
Models:
model.год: value ~ 1 + (1 | Год)
model.ful: value ~ 1 + (1 | Год) + (1 | variable)
          Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)
model.год  3 43.794 44.702 -18.897                        
model.ful  4 44.004 45.215 -18.002 1.7903      1     0.1809
> anova(model.год,model.город)
Data: data.long
Models:
model.год: value ~ 1 + (1 | Год)
model.город: value ~ 1 + (1 | variable)
            Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)    
model.год    3 43.794 44.702 -18.897                            
model.город  3 42.116 43.024 -18.058 1.6782      0  < 2.2e-16 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
>

Сообщение отредактировал p2004r - 2.12.2011 - 14:31

Эскизы прикрепленных изображений