Chap 11
Statistical modeling
March 21, 2017
library(xlsx)## Loading required package: rJava## Loading required package: xlsxjarsa = read.xlsx(file='Table11-8.xls',sheetIndex = 1,header = TRUE)
names(a)[3] = 'y'
names(a)[4] = 'x'
library(ggplot2)
a$dog = as.character(a$dog)
a$condition = as.character(a$condition)
#a
ggplot(data=a,aes(x=condition,y=y,group=dog))+geom_line(aes(colour=dog))
ggplot(data=a,aes(x=condition,y=x,group=dog))+geom_line(aes(colour=dog))
ggplot(data=a,aes(x=x,y=y,group=dog))+geom_point(aes(colour=dog,size=condition))## Warning: Using size for a discrete variable is not advised.
ggplot(data=a,aes(x=x,y=y,group=dog))+geom_point(aes(colour=as.character(condition),shape = dog))
## Make more stuff (geom_boxplot for instance) and analize this full-informative pictures.
#b
mod1 = glm(a$y ~ a$x)
summary(mod1)##
## Call:
## glm(formula = a$y ~ a$x)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -15.8908 -5.3370 0.9632 5.9919 12.9169
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.76808 6.61726 6.161 3.43e-07 ***
## a$x 0.76923 0.09156 8.401 3.42e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 62.57894)
##
## Null deviance: 6795 on 39 degrees of freedom
## Residual deviance: 2378 on 38 degrees of freedom
## AIC: 282.92
##
## Number of Fisher Scoring iterations: 2library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionintercept=NULL
slope=NULL
for(i in 1:5){
b = filter(a,dog==as.character(i))
mod = glm(b$y ~ b$x)
intercept[i] = mod$coefficients[1]
slope[i] = mod$coefficients[2]
}
mean(intercept)## [1] 37.37958mean(slope)## [1] 0.7737112intercept=NULL
slope=NULL
for(i in 1:8){
b = filter(a,condition==as.character(i))
mod = glm(b$y ~ b$x)
intercept[i] = mod$coefficients[1]
slope[i] = mod$coefficients[2]
}
mean(intercept)## [1] 40.25895mean(slope)## [1] 0.7676023#c
library(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse#a$dog = as.integer(a$dog)
mod2 = lme(y ~ x+condition, data=a, random=~1|dog)
summary(mod2)## Linear mixed-effects model fit by REML
## Data: a
## AIC BIC logLik
## 263.8652 279.6391 -120.9326
##
## Random effects:
## Formula: ~1 | dog
## (Intercept) Residual
## StdDev: 0.003109465 8.446683
##
## Fixed effects: y ~ x + condition
## Value Std.Error DF t-value p-value
## (Intercept) 45.39288 8.785049 27 5.167060 0.0000
## x 0.72841 0.107852 27 6.753767 0.0000
## condition2 0.05921 5.358685 27 0.011049 0.9913
## condition3 -3.79530 5.343113 27 -0.710316 0.4836
## condition4 -5.32751 5.472300 27 -0.973542 0.3389
## condition5 -0.81482 5.342214 27 -0.152524 0.8799
## condition6 -0.91219 5.488666 27 -0.166196 0.8692
## condition7 -3.83401 5.372417 27 -0.713646 0.4816
## condition8 0.80069 5.361322 27 0.149345 0.8824
## Correlation:
## (Intr) x cndtn2 cndtn3 cndtn4 cndtn5 cndtn6 cndtn7
## x -0.903
## condition2 -0.232 -0.078
## condition3 -0.321 0.019 0.497
## condition4 -0.493 0.217 0.470 0.492
## condition5 -0.300 -0.005 0.499 0.500 0.487
## condition6 -0.503 0.230 0.467 0.491 0.525 0.486
## condition7 -0.398 0.106 0.487 0.499 0.508 0.497 0.508
## condition8 -0.227 -0.084 0.503 0.497 0.468 0.499 0.466 0.486
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.67845684 -0.74839654 0.05759677 0.59043356 1.57302750
##
## Number of Observations: 40
## Number of Groups: 5mod2 = lme(y ~ x+dog, data=a, random=~1|condition)
summary(mod2)## Linear mixed-effects model fit by REML
## Data: a
## AIC BIC logLik
## 271.514 283.7249 -127.757
##
## Random effects:
## Formula: ~1 | condition
## (Intercept) Residual
## StdDev: 0.0003473208 7.878833
##
## Fixed effects: y ~ x + dog
## Value Std.Error DF t-value p-value
## (Intercept) 44.86257 7.280654 27 6.161888 0.0000
## x 0.62885 0.126412 27 4.974572 0.0000
## dog2 5.12320 5.038551 27 1.016800 0.3183
## dog3 8.17552 4.611402 27 1.772893 0.0875
## dog4 9.07697 5.188606 27 1.749405 0.0916
## dog5 6.96571 4.566030 27 1.525552 0.1388
## Correlation:
## (Intr) x dog2 dog3 dog4
## x -0.924
## dog2 0.365 -0.623
## dog3 0.249 -0.520 0.658
## dog4 0.396 -0.651 0.703 0.663
## dog5 0.234 -0.506 0.653 0.631 0.657
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.7666673 -0.6197425 0.2017892 0.5708534 1.6250549
##
## Number of Observations: 40
## Number of Groups: 8mod2 = lme(y ~ x, data=a, random=~1|condition/dog)
summary(mod2)## Linear mixed-effects model fit by REML
## Data: a
## AIC BIC logLik
## 287.6305 295.8184 -138.8152
##
## Random effects:
## Formula: ~1 | condition
## (Intercept)
## StdDev: 0.0002590303
##
## Formula: ~1 | dog %in% condition
## (Intercept) Residual
## StdDev: 7.91057 0.0427298
##
## Fixed effects: y ~ x
## Value Std.Error DF t-value p-value
## (Intercept) 40.76808 6.617261 31 6.160870 0
## x 0.76923 0.091559 31 8.401375 0
## Correlation:
## (Intr)
## x -0.982
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.0108504785 -0.0036441991 0.0006576696 0.0040913580 0.0088198513
##
## Number of Observations: 40
## Number of Groups:
## condition dog %in% condition
## 8 40mod2 = lme(y ~ x, data=a, random=~1|dog/condition)
summary(mod2)## Linear mixed-effects model fit by REML
## Data: a
## AIC BIC logLik
## 287.6305 295.8184 -138.8152
##
## Random effects:
## Formula: ~1 | dog
## (Intercept)
## StdDev: 0.0008065215
##
## Formula: ~1 | condition %in% dog
## (Intercept) Residual
## StdDev: 7.910363 0.07147414
##
## Fixed effects: y ~ x
## Value Std.Error DF t-value p-value
## (Intercept) 40.76808 6.617261 34 6.160870 0
## x 0.76923 0.091559 34 8.401375 0
## Correlation:
## (Intr)
## x -0.982
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.018149595 -0.006095652 0.001100084 0.006843614 0.014752965
##
## Number of Observations: 40
## Number of Groups:
## dog condition %in% dog
## 5 40library(geepack)
mod2 = geeglm(y ~ x+condition, data=a, id=dog, wave=x)
summary(mod2)##
## Call:
## geeglm(formula = y ~ x + condition, data = a, id = dog, waves = x)
##
## Coefficients:
## Estimate Std.err Wald Pr(>|W|)
## (Intercept) 45.39287 6.77520 44.888 2.09e-11 ***
## x 0.72841 0.08703 70.057 < 2e-16 ***
## condition2 0.05921 2.22420 0.001 0.9788
## condition3 -3.79530 3.97435 0.912 0.3396
## condition4 -5.32751 2.95438 3.252 0.0713 .
## condition5 -0.81482 4.29508 0.036 0.8495
## condition6 -0.91219 3.66640 0.062 0.8035
## condition7 -3.83401 4.37933 0.766 0.3813
## condition8 0.80069 3.92078 0.042 0.8382
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Estimated Scale Parameters:
## Estimate Std.err
## (Intercept) 55.29 10.2
##
## Correlation: Structure = independenceNumber of clusters: 5 Maximum cluster size: 8mod2 = geeglm(y ~ x+dog, data=a, id=condition, wave=x)
summary(mod2)##
## Call:
## geeglm(formula = y ~ x + dog, data = a, id = condition, waves = x)
##
## Coefficients:
## Estimate Std.err Wald Pr(>|W|)
## (Intercept) 44.8626 4.6430 93.36 < 2e-16 ***
## x 0.6288 0.0818 59.04 1.5e-14 ***
## dog2 5.1232 4.7041 1.19 0.276
## dog3 8.1755 3.5073 5.43 0.020 *
## dog4 9.0770 4.0122 5.12 0.024 *
## dog5 6.9657 3.5274 3.90 0.048 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Estimated Scale Parameters:
## Estimate Std.err
## (Intercept) 52.8 8.86
##
## Correlation: Structure = independenceNumber of clusters: 40 Maximum cluster size: 1```