reshape_plyr.R

etienne — Jun 6, 2013, 11:09 PM

# reshape and plyr data manipulation tutorial
# Etienne Laliberte
# June 5, 2013


# load libraries
library(reshape2)
library(plyr)

###############################################################
# reshape


# the smiths data
# convert time to a factor
smiths$time <- as.factor(smiths$time)


### melting the data
melt(smiths)
Using subject, time as id variables
     subject time variable value
1 John Smith    1      age 33.00
2 Mary Smith    1      age    NA
3 John Smith    1   weight 90.00
4 Mary Smith    1   weight    NA
5 John Smith    1   height  1.87
6 Mary Smith    1   height  1.54
# the following does the same thing
# just different ways of manually specifying
# which variables are identifiers and which are measured
melt(smiths, id = c("subject", "time"), measured = c("age", "weight",
                                                     "height"))
     subject time variable value
1 John Smith    1      age 33.00
2 Mary Smith    1      age    NA
3 John Smith    1   weight 90.00
4 Mary Smith    1   weight    NA
5 John Smith    1   height  1.87
6 Mary Smith    1   height  1.54
melt(smiths, id = c("subject", "time"))
     subject time variable value
1 John Smith    1      age 33.00
2 Mary Smith    1      age    NA
3 John Smith    1   weight 90.00
4 Mary Smith    1   weight    NA
5 John Smith    1   height  1.87
6 Mary Smith    1   height  1.54
melt(smiths, id = 1:2)
     subject time variable value
1 John Smith    1      age 33.00
2 Mary Smith    1      age    NA
3 John Smith    1   weight 90.00
4 Mary Smith    1   weight    NA
5 John Smith    1   height  1.87
6 Mary Smith    1   height  1.54
melt(smiths, measured = c("age", "weight", "height"))
Using subject, time as id variables
     subject time variable value
1 John Smith    1      age 33.00
2 Mary Smith    1      age    NA
3 John Smith    1   weight 90.00
4 Mary Smith    1   weight    NA
5 John Smith    1   height  1.87
6 Mary Smith    1   height  1.54


### structural missings
melt(smiths, na.rm = TRUE)
Using subject, time as id variables
     subject time variable value
1 John Smith    1      age 33.00
3 John Smith    1   weight 90.00
5 John Smith    1   height  1.87
6 Mary Smith    1   height  1.54


# casting the data
smithsm <- melt(smiths, na.rm = TRUE)
Using subject, time as id variables
dcast(smithsm, time + subject ~ variable)
  time    subject age weight height
1    1 John Smith  33     90   1.87
2    1 Mary Smith  NA     NA   1.54
smiths # same as original data
     subject time age weight height
1 John Smith    1  33     90   1.87
2 Mary Smith    1  NA     NA   1.54


# casting formulas
dcast(smithsm, subject ~ .)
Aggregation function missing: defaulting to length
     subject NA
1 John Smith  3
2 Mary Smith  1
dcast(smithsm, ... ~ variable)
     subject time age weight height
1 John Smith    1  33     90   1.87
2 Mary Smith    1  NA     NA   1.54
dcast(smithsm, time + subject ~ variable)
  time    subject age weight height
1    1 John Smith  33     90   1.87
2    1 Mary Smith  NA     NA   1.54
dcast(smithsm, ... ~ subject)
  time variable John Smith Mary Smith
1    1      age      33.00         NA
2    1   weight      90.00         NA
3    1   height       1.87       1.54
dcast(smithsm, ... ~ time)
     subject variable     1
1 John Smith      age 33.00
2 John Smith   weight 90.00
3 John Smith   height  1.87
4 Mary Smith   height  1.54


# root layer data
rootdata <- read.csv("root layer data.csv")
# melt data frame
# add pot factor
rootdata$ID <- as.factor(rootdata$ID)
(rootdata.melt <- melt(rootdata))[1:3,]
Using ID, P, W as id variables
  ID   P   W variable value
1  1 hiP hiW   r10.h1 1.979
2  2 hiP hiW   r10.h1 2.921
3  3 hiP hiW   r10.h1 1.797

# create new variables depth and harvest
rootdata.melt$layer <- as.factor(substr(rootdata.melt$variable, start = 1, stop = 3) )
rootdata.melt$harvest <- as.factor(substr(rootdata.melt$variable, start = 5, stop = 6) )
rootdata.melt[1:3,]
  ID   P   W variable value layer harvest
1  1 hiP hiW   r10.h1 1.979   r10      h1
2  2 hiP hiW   r10.h1 2.921   r10      h1
3  3 hiP hiW   r10.h1 1.797   r10      h1


# aggregating the data
# how many replicates of each treatment?
table(rootdata$P,rootdata$W)

       hiW loW MedW
  hiP    6   3    6
  loP    3   6    6
  medP   6   6    3
dcast(rootdata.melt, P + W + layer + harvest ~ .)
Aggregation function missing: defaulting to length
      P    W layer harvest NA
1   hiP  hiW   r10      h1  6
2   hiP  hiW   r10      h2  6
3   hiP  hiW   r20      h1  6
4   hiP  hiW   r20      h2  6
5   hiP  hiW   r30      h1  6
6   hiP  hiW   r30      h2  6
7   hiP  hiW   r40      h1  6
8   hiP  hiW   r40      h2  6
9   hiP  hiW   r50      h1  6
10  hiP  hiW   r50      h2  6
11  hiP  loW   r10      h1  3
12  hiP  loW   r10      h2  3
13  hiP  loW   r20      h1  3
14  hiP  loW   r20      h2  3
15  hiP  loW   r30      h1  3
16  hiP  loW   r30      h2  3
17  hiP  loW   r40      h1  3
18  hiP  loW   r40      h2  3
19  hiP  loW   r50      h1  3
20  hiP  loW   r50      h2  3
21  hiP MedW   r10      h1  6
22  hiP MedW   r10      h2  6
23  hiP MedW   r20      h1  6
24  hiP MedW   r20      h2  6
25  hiP MedW   r30      h1  6
26  hiP MedW   r30      h2  6
27  hiP MedW   r40      h1  6
28  hiP MedW   r40      h2  6
29  hiP MedW   r50      h1  6
30  hiP MedW   r50      h2  6
31  loP  hiW   r10      h1  3
32  loP  hiW   r10      h2  3
33  loP  hiW   r20      h1  3
34  loP  hiW   r20      h2  3
35  loP  hiW   r30      h1  3
36  loP  hiW   r30      h2  3
37  loP  hiW   r40      h1  3
38  loP  hiW   r40      h2  3
39  loP  hiW   r50      h1  3
40  loP  hiW   r50      h2  3
41  loP  loW   r10      h1  6
42  loP  loW   r10      h2  6
43  loP  loW   r20      h1  6
44  loP  loW   r20      h2  6
45  loP  loW   r30      h1  6
46  loP  loW   r30      h2  6
47  loP  loW   r40      h1  6
48  loP  loW   r40      h2  6
49  loP  loW   r50      h1  6
50  loP  loW   r50      h2  6
51  loP MedW   r10      h1  6
52  loP MedW   r10      h2  6
53  loP MedW   r20      h1  6
54  loP MedW   r20      h2  6
55  loP MedW   r30      h1  6
56  loP MedW   r30      h2  6
57  loP MedW   r40      h1  6
58  loP MedW   r40      h2  6
59  loP MedW   r50      h1  6
60  loP MedW   r50      h2  6
61 medP  hiW   r10      h1  6
62 medP  hiW   r10      h2  6
63 medP  hiW   r20      h1  6
64 medP  hiW   r20      h2  6
65 medP  hiW   r30      h1  6
66 medP  hiW   r30      h2  6
67 medP  hiW   r40      h1  6
68 medP  hiW   r40      h2  6
69 medP  hiW   r50      h1  6
70 medP  hiW   r50      h2  6
71 medP  loW   r10      h1  6
72 medP  loW   r10      h2  6
73 medP  loW   r20      h1  6
74 medP  loW   r20      h2  6
75 medP  loW   r30      h1  6
76 medP  loW   r30      h2  6
77 medP  loW   r40      h1  6
78 medP  loW   r40      h2  6
79 medP  loW   r50      h1  6
80 medP  loW   r50      h2  6
81 medP MedW   r10      h1  3
82 medP MedW   r10      h2  3
83 medP MedW   r20      h1  3
84 medP MedW   r20      h2  3
85 medP MedW   r30      h1  3
86 medP MedW   r30      h2  3
87 medP MedW   r40      h1  3
88 medP MedW   r40      h2  3
89 medP MedW   r50      h1  3
90 medP MedW   r50      h2  3

# calculate the medians of each measured variable for each treatement
# approach you saw this morning
aggregate(rootdata[, 4:13], by = list(rootdata$P, rootdata$W), median) 
  Group.1 Group.2 r10.h1 r20.h1 r30.h1 r40.h1 r50.h1 r10.h2 r50.h2 r20.h2
1     hiP     hiW  1.912 2.0680 0.0690      0      0  3.365  2.546 0.3975
2     loP     hiW  2.314 0.3970 0.0850      0      0  3.414  0.920 0.0820
3    medP     hiW  2.230 0.9220 0.1115      0      0  4.826  1.784 0.2235
4     hiP     loW  1.464 1.2550 0.0120      0      0  1.009  3.697 0.3450
5     loP     loW  1.565 0.3940 0.0705      0      0  2.071  1.332 0.1570
6    medP     loW  1.966 0.8340 0.0780      0      0  2.166  2.533 0.2850
7     hiP    MedW  1.837 1.7250 0.0435      0      0  2.435  3.964 0.5055
8     loP    MedW  1.580 0.3895 0.0810      0      0  3.464  1.447 0.1865
9    medP    MedW  1.444 0.8480 0.0990      0      0  2.490  2.407 0.4200
  r30.h2 r40.h2
1 0.2420 0.1520
2 0.0870 0.0300
3 0.1840 0.0595
4 0.1990 0.1180
5 0.1270 0.0645
6 0.2720 0.1270
7 0.3545 0.2435
8 0.1105 0.0585
9 0.2740 0.0810
# with dcast
dcast(rootdata.melt, P + W ~  harvest + layer, median)
     P    W h1_r10 h1_r20 h1_r30 h1_r40 h1_r50 h2_r10 h2_r20 h2_r30 h2_r40
1  hiP  hiW  1.912 2.0680 0.0690      0      0  3.365 0.3975 0.2420 0.1520
2  hiP  loW  1.464 1.2550 0.0120      0      0  1.009 0.3450 0.1990 0.1180
3  hiP MedW  1.837 1.7250 0.0435      0      0  2.435 0.5055 0.3545 0.2435
4  loP  hiW  2.314 0.3970 0.0850      0      0  3.414 0.0820 0.0870 0.0300
5  loP  loW  1.565 0.3940 0.0705      0      0  2.071 0.1570 0.1270 0.0645
6  loP MedW  1.580 0.3895 0.0810      0      0  3.464 0.1865 0.1105 0.0585
7 medP  hiW  2.230 0.9220 0.1115      0      0  4.826 0.2235 0.1840 0.0595
8 medP  loW  1.966 0.8340 0.0780      0      0  2.166 0.2850 0.2720 0.1270
9 medP MedW  1.444 0.8480 0.0990      0      0  2.490 0.4200 0.2740 0.0810
  h2_r50
1  2.546
2  3.697
3  3.964
4  0.920
5  1.332
6  1.447
7  1.784
8  2.533
9  2.407

# calculate total root length per pot
# approch you saw this morning
rowSums(rootdata[,4:8])
 [1] 4.218 4.984 3.288 4.183 4.095 3.646 3.212 3.101 2.517 3.853 3.849
[12] 5.204 3.576 4.023 3.088 2.994 3.078 2.925 3.132 4.211 3.079 2.272
[23] 2.342 3.125 2.307 3.025 2.219 4.218 3.122 3.356 2.539 1.879 2.208
[34] 1.766 1.259 2.278 2.633 3.060 2.888 1.784 2.522 1.805 2.065 1.919
[45] 2.355
rowSums(rootdata[,9:13])
 [1]  6.246  6.630  5.692  7.401  8.112 10.264  5.053  7.230  5.041  8.338
[11]  8.816  7.214  8.389  7.642  6.841  5.304  7.063  4.391  9.137  7.148
[21]  8.729  6.239  5.559  6.964  5.743  5.894  3.837  5.799  8.578  6.368
[31]  3.954  3.313  5.470  6.408  3.699  4.809  9.381  3.749  4.746  8.797
[41]  7.491  4.855  3.707  3.710  2.735
# with dcast
dcast(rootdata.melt, P + W + ID ~ harvest, sum)
      P    W ID    h1     h2
1   hiP  hiW  1 4.218  6.246
2   hiP  hiW  2 4.984  6.630
3   hiP  hiW  3 3.288  5.692
4   hiP  hiW 10 3.853  8.338
5   hiP  hiW 11 3.849  8.816
6   hiP  hiW 12 5.204  7.214
7   hiP  loW  7 3.212  5.053
8   hiP  loW  8 3.101  7.230
9   hiP  loW  9 2.517  5.041
10  hiP MedW  4 4.183  7.401
11  hiP MedW  5 4.095  8.112
12  hiP MedW  6 3.646 10.264
13  hiP MedW 13 3.576  8.389
14  hiP MedW 14 4.023  7.642
15  hiP MedW 15 3.088  6.841
16  loP  hiW 37 2.633  9.381
17  loP  hiW 38 3.060  3.749
18  loP  hiW 39 2.888  4.746
19  loP  loW 34 1.766  6.408
20  loP  loW 35 1.259  3.699
21  loP  loW 36 2.278  4.809
22  loP  loW 43 2.065  3.707
23  loP  loW 44 1.919  3.710
24  loP  loW 45 2.355  2.735
25  loP MedW 31 2.539  3.954
26  loP MedW 32 1.879  3.313
27  loP MedW 33 2.208  5.470
28  loP MedW 40 1.784  8.797
29  loP MedW 41 2.522  7.491
30  loP MedW 42 1.805  4.855
31 medP  hiW 19 3.132  9.137
32 medP  hiW 20 4.211  7.148
33 medP  hiW 21 3.079  8.729
34 medP  hiW 28 4.218  5.799
35 medP  hiW 29 3.122  8.578
36 medP  hiW 30 3.356  6.368
37 medP  loW 16 2.994  5.304
38 medP  loW 17 3.078  7.063
39 medP  loW 18 2.925  4.391
40 medP  loW 25 2.307  5.743
41 medP  loW 26 3.025  5.894
42 medP  loW 27 2.219  3.837
43 medP MedW 22 2.272  6.239
44 medP MedW 23 2.342  5.559
45 medP MedW 24 3.125  6.964


# exercices
# 1) mean root length per layer for each W/P combination across both harvests
# 2) sd of different layers for each water treatment per harvest


# margins
dcast(rootdata.melt, W ~ P, mean, margins = TRUE)
      W    hiP    loP   medP  (all)
1   hiW 1.1389 0.8819 1.1146 1.0778
2   loW 0.8718 0.6118 0.8130 0.7443
3  MedW 1.1877 0.7770 0.8834 0.9625
4 (all) 1.1050 0.7319 0.9477 0.9282
dcast(rootdata.melt, W ~ P, mean, margins = 'P')
     W    hiP    loP   medP  (all)
1  hiW 1.1389 0.8819 1.1146 1.0778
2  loW 0.8718 0.6118 0.8130 0.7443
3 MedW 1.1877 0.7770 0.8834 0.9625
dcast(rootdata.melt, W ~ P, mean, margins = c('W', 'P') ) # same as margins = T
      W    hiP    loP   medP  (all)
1   hiW 1.1389 0.8819 1.1146 1.0778
2   loW 0.8718 0.6118 0.8130 0.7443
3  MedW 1.1877 0.7770 0.8834 0.9625
4 (all) 1.1050 0.7319 0.9477 0.9282



###################################################################
#### plyr


# group-wise ranking
rootdata.melt[1:20,]
   ID    P    W variable value layer harvest
1   1  hiP  hiW   r10.h1 1.979   r10      h1
2   2  hiP  hiW   r10.h1 2.921   r10      h1
3   3  hiP  hiW   r10.h1 1.797   r10      h1
4   4  hiP MedW   r10.h1 1.950   r10      h1
5   5  hiP MedW   r10.h1 2.608   r10      h1
6   6  hiP MedW   r10.h1 1.901   r10      h1
7   7  hiP  loW   r10.h1 1.464   r10      h1
8   8  hiP  loW   r10.h1 1.834   r10      h1
9   9  hiP  loW   r10.h1 1.172   r10      h1
10 10  hiP  hiW   r10.h1 1.270   r10      h1
11 11  hiP  hiW   r10.h1 1.844   r10      h1
12 12  hiP  hiW   r10.h1 2.901   r10      h1
13 13  hiP MedW   r10.h1 1.772   r10      h1
14 14  hiP MedW   r10.h1 1.650   r10      h1
15 15  hiP MedW   r10.h1 1.327   r10      h1
16 16 medP  loW   r10.h1 2.067   r10      h1
17 17 medP  loW   r10.h1 2.016   r10      h1
18 18 medP  loW   r10.h1 2.269   r10      h1
19 19 medP  hiW   r10.h1 1.921   r10      h1
20 20 medP  hiW   r10.h1 3.302   r10      h1


# 1. extract a piece of the data
rootsub <- subset(rootdata.melt, P == 'hiP' & W == 'hiW' & layer == 'r10' & harvest == 'h1')
rootsub
   ID   P   W variable value layer harvest
1   1 hiP hiW   r10.h1 1.979   r10      h1
2   2 hiP hiW   r10.h1 2.921   r10      h1
3   3 hiP hiW   r10.h1 1.797   r10      h1
10 10 hiP hiW   r10.h1 1.270   r10      h1
11 11 hiP hiW   r10.h1 1.844   r10      h1
12 12 hiP hiW   r10.h1 2.901   r10      h1

# 2. solve the problem by hand
rank(rootsub$value)
[1] 4 6 2 1 3 5
rootsub$ranks <- rank(rootsub$value)
rootsub
   ID   P   W variable value layer harvest ranks
1   1 hiP hiW   r10.h1 1.979   r10      h1     4
2   2 hiP hiW   r10.h1 2.921   r10      h1     6
3   3 hiP hiW   r10.h1 1.797   r10      h1     2
10 10 hiP hiW   r10.h1 1.270   r10      h1     1
11 11 hiP hiW   r10.h1 1.844   r10      h1     3
12 12 hiP hiW   r10.h1 2.901   r10      h1     5


# 3. write a function that solves the problem
myfun <- function(x) {
  x$ranks <- rank(x$value)
  return(x)
}


# 4. use plyr to perform operation on each group
myranks <- ddply(rootdata.melt, .(P, W, layer, harvest), myfun)
myranks[1:12,]
   ID   P   W variable value layer harvest ranks
1   1 hiP hiW   r10.h1 1.979   r10      h1     4
2   2 hiP hiW   r10.h1 2.921   r10      h1     6
3   3 hiP hiW   r10.h1 1.797   r10      h1     2
4  10 hiP hiW   r10.h1 1.270   r10      h1     1
5  11 hiP hiW   r10.h1 1.844   r10      h1     3
6  12 hiP hiW   r10.h1 2.901   r10      h1     5
7   1 hiP hiW   r10.h2 3.391   r10      h2     4
8   2 hiP hiW   r10.h2 2.408   r10      h2     1
9   3 hiP hiW   r10.h2 2.631   r10      h2     2
10 10 hiP hiW   r10.h2 3.340   r10      h2     3
11 11 hiP hiW   r10.h2 5.163   r10      h2     6
12 12 hiP hiW   r10.h2 3.836   r10      h2     5
# this is the same
myranks <- ddply(rootdata.melt, c('P', 'W', 'layer', 'harvest'), myfun) # characters
myranks <- ddply(rootdata.melt, ~ P + W + layer + harvest, myfun) # formula


# column-wise functions
ddply(rootdata.melt, .(P, W, layer, harvest), numcolwise(rank ))
       P    W layer harvest value
1    hiP  hiW   r10      h1   4.0
2    hiP  hiW   r10      h1   6.0
3    hiP  hiW   r10      h1   2.0
4    hiP  hiW   r10      h1   1.0
5    hiP  hiW   r10      h1   3.0
6    hiP  hiW   r10      h1   5.0
7    hiP  hiW   r10      h2   4.0
8    hiP  hiW   r10      h2   1.0
9    hiP  hiW   r10      h2   2.0
10   hiP  hiW   r10      h2   3.0
11   hiP  hiW   r10      h2   6.0
12   hiP  hiW   r10      h2   5.0
13   hiP  hiW   r20      h1   4.0
14   hiP  hiW   r20      h1   3.0
15   hiP  hiW   r20      h1   1.0
16   hiP  hiW   r20      h1   6.0
17   hiP  hiW   r20      h1   2.0
18   hiP  hiW   r20      h1   5.0
19   hiP  hiW   r20      h2   2.5
20   hiP  hiW   r20      h2   1.0
21   hiP  hiW   r20      h2   2.5
22   hiP  hiW   r20      h2   5.0
23   hiP  hiW   r20      h2   6.0
24   hiP  hiW   r20      h2   4.0
25   hiP  hiW   r30      h1   5.0
26   hiP  hiW   r30      h1   2.0
27   hiP  hiW   r30      h1   1.0
28   hiP  hiW   r30      h1   4.0
29   hiP  hiW   r30      h1   6.0
30   hiP  hiW   r30      h1   3.0
31   hiP  hiW   r30      h2   4.0
32   hiP  hiW   r30      h2   1.0
33   hiP  hiW   r30      h2   3.0
34   hiP  hiW   r30      h2   5.0
35   hiP  hiW   r30      h2   2.0
36   hiP  hiW   r30      h2   6.0
37   hiP  hiW   r40      h1   3.5
38   hiP  hiW   r40      h1   3.5
39   hiP  hiW   r40      h1   3.5
40   hiP  hiW   r40      h1   3.5
41   hiP  hiW   r40      h1   3.5
42   hiP  hiW   r40      h1   3.5
43   hiP  hiW   r40      h2   3.0
44   hiP  hiW   r40      h2   4.0
45   hiP  hiW   r40      h2   5.0
46   hiP  hiW   r40      h2   6.0
47   hiP  hiW   r40      h2   1.0
48   hiP  hiW   r40      h2   2.0
49   hiP  hiW   r50      h1   3.5
50   hiP  hiW   r50      h1   3.5
51   hiP  hiW   r50      h1   3.5
52   hiP  hiW   r50      h1   3.5
53   hiP  hiW   r50      h1   3.5
54   hiP  hiW   r50      h1   3.5
55   hiP  hiW   r50      h2   1.0
56   hiP  hiW   r50      h2   5.0
57   hiP  hiW   r50      h2   2.0
58   hiP  hiW   r50      h2   6.0
59   hiP  hiW   r50      h2   4.0
60   hiP  hiW   r50      h2   3.0
61   hiP  loW   r10      h1   2.0
62   hiP  loW   r10      h1   3.0
63   hiP  loW   r10      h1   1.0
64   hiP  loW   r10      h2   2.0
65   hiP  loW   r10      h2   3.0
66   hiP  loW   r10      h2   1.0
67   hiP  loW   r20      h1   3.0
68   hiP  loW   r20      h1   2.0
69   hiP  loW   r20      h1   1.0
70   hiP  loW   r20      h2   3.0
71   hiP  loW   r20      h2   1.0
72   hiP  loW   r20      h2   2.0
73   hiP  loW   r30      h1   1.0
74   hiP  loW   r30      h1   2.0
75   hiP  loW   r30      h1   3.0
76   hiP  loW   r30      h2   3.0
77   hiP  loW   r30      h2   2.0
78   hiP  loW   r30      h2   1.0
79   hiP  loW   r40      h1   2.0
80   hiP  loW   r40      h1   2.0
81   hiP  loW   r40      h1   2.0
82   hiP  loW   r40      h2   1.0
83   hiP  loW   r40      h2   3.0
84   hiP  loW   r40      h2   2.0
85   hiP  loW   r50      h1   2.0
86   hiP  loW   r50      h1   2.0
87   hiP  loW   r50      h1   2.0
88   hiP  loW   r50      h2   1.0
89   hiP  loW   r50      h2   3.0
90   hiP  loW   r50      h2   2.0
91   hiP MedW   r10      h1   5.0
92   hiP MedW   r10      h1   6.0
93   hiP MedW   r10      h1   4.0
94   hiP MedW   r10      h1   3.0
95   hiP MedW   r10      h1   2.0
96   hiP MedW   r10      h1   1.0
97   hiP MedW   r10      h2   4.0
98   hiP MedW   r10      h2   6.0
99   hiP MedW   r10      h2   3.0
100  hiP MedW   r10      h2   1.0
101  hiP MedW   r10      h2   5.0
102  hiP MedW   r10      h2   2.0
103  hiP MedW   r20      h1   5.0
104  hiP MedW   r20      h1   1.0
105  hiP MedW   r20      h1   3.0
106  hiP MedW   r20      h1   2.0
107  hiP MedW   r20      h1   6.0
108  hiP MedW   r20      h1   4.0
109  hiP MedW   r20      h2   4.0
110  hiP MedW   r20      h2   2.0
111  hiP MedW   r20      h2   6.0
112  hiP MedW   r20      h2   5.0
113  hiP MedW   r20      h2   1.0
114  hiP MedW   r20      h2   3.0
115  hiP MedW   r30      h1   4.0
116  hiP MedW   r30      h1   1.0
117  hiP MedW   r30      h1   3.0
118  hiP MedW   r30      h1   6.0
119  hiP MedW   r30      h1   5.0
120  hiP MedW   r30      h1   2.0
121  hiP MedW   r30      h2   1.0
122  hiP MedW   r30      h2   5.0
123  hiP MedW   r30      h2   6.0
124  hiP MedW   r30      h2   3.0
125  hiP MedW   r30      h2   4.0
126  hiP MedW   r30      h2   2.0
127  hiP MedW   r40      h1   3.5
128  hiP MedW   r40      h1   3.5
129  hiP MedW   r40      h1   3.5
130  hiP MedW   r40      h1   3.5
131  hiP MedW   r40      h1   3.5
132  hiP MedW   r40      h1   3.5
133  hiP MedW   r40      h2   5.0
134  hiP MedW   r40      h2   6.0
135  hiP MedW   r40      h2   4.0
136  hiP MedW   r40      h2   3.0
137  hiP MedW   r40      h2   1.0
138  hiP MedW   r40      h2   2.0
139  hiP MedW   r50      h1   3.5
140  hiP MedW   r50      h1   3.5
141  hiP MedW   r50      h1   3.5
142  hiP MedW   r50      h1   3.5
143  hiP MedW   r50      h1   3.5
144  hiP MedW   r50      h1   3.5
145  hiP MedW   r50      h2   3.0
146  hiP MedW   r50      h2   1.0
147  hiP MedW   r50      h2   6.0
148  hiP MedW   r50      h2   5.0
149  hiP MedW   r50      h2   4.0
150  hiP MedW   r50      h2   2.0
151  loP  hiW   r10      h1   1.0
152  loP  hiW   r10      h1   3.0
153  loP  hiW   r10      h1   2.0
154  loP  hiW   r10      h2   3.0
155  loP  hiW   r10      h2   1.0
156  loP  hiW   r10      h2   2.0
157  loP  hiW   r20      h1   1.0
158  loP  hiW   r20      h1   2.0
159  loP  hiW   r20      h1   3.0
160  loP  hiW   r20      h2   3.0
161  loP  hiW   r20      h2   2.0
162  loP  hiW   r20      h2   1.0
163  loP  hiW   r30      h1   1.0
164  loP  hiW   r30      h1   3.0
165  loP  hiW   r30      h1   2.0
166  loP  hiW   r30      h2   1.0
167  loP  hiW   r30      h2   3.0
168  loP  hiW   r30      h2   2.0
169  loP  hiW   r40      h1   2.0
170  loP  hiW   r40      h1   2.0
171  loP  hiW   r40      h1   2.0
172  loP  hiW   r40      h2   1.0
173  loP  hiW   r40      h2   2.0
174  loP  hiW   r40      h2   3.0
175  loP  hiW   r50      h1   2.0
176  loP  hiW   r50      h1   2.0
177  loP  hiW   r50      h1   2.0
178  loP  hiW   r50      h2   1.0
179  loP  hiW   r50      h2   2.0
180  loP  hiW   r50      h2   3.0
181  loP  loW   r10      h1   3.0
182  loP  loW   r10      h1   1.0
183  loP  loW   r10      h1   4.0
184  loP  loW   r10      h1   2.0
185  loP  loW   r10      h1   5.0
186  loP  loW   r10      h1   6.0
187  loP  loW   r10      h2   5.0
188  loP  loW   r10      h2   4.0
189  loP  loW   r10      h2   6.0
190  loP  loW   r10      h2   3.0
191  loP  loW   r10      h2   1.5
192  loP  loW   r10      h2   1.5
193  loP  loW   r20      h1   2.0
194  loP  loW   r20      h1   3.5
195  loP  loW   r20      h1   5.0
196  loP  loW   r20      h1   6.0
197  loP  loW   r20      h1   1.0
198  loP  loW   r20      h1   3.5
199  loP  loW   r20      h2   3.0
200  loP  loW   r20      h2   4.0
201  loP  loW   r20      h2   6.0
202  loP  loW   r20      h2   5.0
203  loP  loW   r20      h2   1.0
204  loP  loW   r20      h2   2.0
205  loP  loW   r30      h1   2.0
206  loP  loW   r30      h1   1.0
207  loP  loW   r30      h1   6.0
208  loP  loW   r30      h1   5.0
209  loP  loW   r30      h1   4.0
210  loP  loW   r30      h1   3.0
211  loP  loW   r30      h2   3.0
212  loP  loW   r30      h2   5.0
213  loP  loW   r30      h2   1.0
214  loP  loW   r30      h2   4.0
215  loP  loW   r30      h2   6.0
216  loP  loW   r30      h2   2.0
217  loP  loW   r40      h1   3.5
218  loP  loW   r40      h1   3.5
219  loP  loW   r40      h1   3.5
220  loP  loW   r40      h1   3.5
221  loP  loW   r40      h1   3.5
222  loP  loW   r40      h1   3.5
223  loP  loW   r40      h2   2.0
224  loP  loW   r40      h2   6.0
225  loP  loW   r40      h2   5.0
226  loP  loW   r40      h2   3.0
227  loP  loW   r40      h2   4.0
228  loP  loW   r40      h2   1.0
229  loP  loW   r50      h1   3.5
230  loP  loW   r50      h1   3.5
231  loP  loW   r50      h1   3.5
232  loP  loW   r50      h1   3.5
233  loP  loW   r50      h1   3.5
234  loP  loW   r50      h1   3.5
235  loP  loW   r50      h2   6.0
236  loP  loW   r50      h2   2.0
237  loP  loW   r50      h2   3.0
238  loP  loW   r50      h2   4.0
239  loP  loW   r50      h2   5.0
240  loP  loW   r50      h2   1.0
241  loP MedW   r10      h1   6.0
242  loP MedW   r10      h1   3.0
243  loP MedW   r10      h1   4.0
244  loP MedW   r10      h1   1.0
245  loP MedW   r10      h1   5.0
246  loP MedW   r10      h1   2.0
247  loP MedW   r10      h2   2.0
248  loP MedW   r10      h2   1.0
249  loP MedW   r10      h2   4.0
250  loP MedW   r10      h2   6.0
251  loP MedW   r10      h2   5.0
252  loP MedW   r10      h2   3.0
253  loP MedW   r20      h1   1.0
254  loP MedW   r20      h1   5.0
255  loP MedW   r20      h1   4.0
256  loP MedW   r20      h1   6.0
257  loP MedW   r20      h1   2.0
258  loP MedW   r20      h1   3.0
259  loP MedW   r20      h2   1.0
260  loP MedW   r20      h2   4.0
261  loP MedW   r20      h2   2.0
262  loP MedW   r20      h2   6.0
263  loP MedW   r20      h2   3.0
264  loP MedW   r20      h2   5.0
265  loP MedW   r30      h1   4.0
266  loP MedW   r30      h1   1.0
267  loP MedW   r30      h1   6.0
268  loP MedW   r30      h1   3.0
269  loP MedW   r30      h1   5.0
270  loP MedW   r30      h1   2.0
271  loP MedW   r30      h2   1.0
272  loP MedW   r30      h2   3.0
273  loP MedW   r30      h2   2.0
274  loP MedW   r30      h2   4.0
275  loP MedW   r30      h2   5.0
276  loP MedW   r30      h2   6.0
277  loP MedW   r40      h1   3.5
278  loP MedW   r40      h1   3.5
279  loP MedW   r40      h1   3.5
280  loP MedW   r40      h1   3.5
281  loP MedW   r40      h1   3.5
282  loP MedW   r40      h1   3.5
283  loP MedW   r40      h2   5.0
284  loP MedW   r40      h2   2.0
285  loP MedW   r40      h2   1.0
286  loP MedW   r40      h2   6.0
287  loP MedW   r40      h2   3.0
288  loP MedW   r40      h2   4.0
289  loP MedW   r50      h1   3.5
290  loP MedW   r50      h1   3.5
291  loP MedW   r50      h1   3.5
292  loP MedW   r50      h1   3.5
293  loP MedW   r50      h1   3.5
294  loP MedW   r50      h1   3.5
295  loP MedW   r50      h2   3.0
296  loP MedW   r50      h2   2.0
297  loP MedW   r50      h2   1.0
298  loP MedW   r50      h2   6.0
299  loP MedW   r50      h2   4.0
300  loP MedW   r50      h2   5.0
301 medP  hiW   r10      h1   1.0
302 medP  hiW   r10      h1   6.0
303 medP  hiW   r10      h1   2.0
304 medP  hiW   r10      h1   5.0
305 medP  hiW   r10      h1   3.0
306 medP  hiW   r10      h1   4.0
307 medP  hiW   r10      h2   3.0
308 medP  hiW   r10      h2   4.0
309 medP  hiW   r10      h2   5.0
310 medP  hiW   r10      h2   2.0
311 medP  hiW   r10      h2   6.0
312 medP  hiW   r10      h2   1.0
313 medP  hiW   r20      h1   6.0
314 medP  hiW   r20      h1   1.0
315 medP  hiW   r20      h1   2.0
316 medP  hiW   r20      h1   4.0
317 medP  hiW   r20      h1   3.0
318 medP  hiW   r20      h1   5.0
319 medP  hiW   r20      h2   3.0
320 medP  hiW   r20      h2   6.0
321 medP  hiW   r20      h2   1.0
322 medP  hiW   r20      h2   2.0
323 medP  hiW   r20      h2   5.0
324 medP  hiW   r20      h2   4.0
325 medP  hiW   r30      h1   3.0
326 medP  hiW   r30      h1   2.0
327 medP  hiW   r30      h1   5.0
328 medP  hiW   r30      h1   6.0
329 medP  hiW   r30      h1   1.0
330 medP  hiW   r30      h1   4.0
331 medP  hiW   r30      h2   6.0
332 medP  hiW   r30      h2   2.0
333 medP  hiW   r30      h2   5.0
334 medP  hiW   r30      h2   1.0
335 medP  hiW   r30      h2   3.0
336 medP  hiW   r30      h2   4.0
337 medP  hiW   r40      h1   3.5
338 medP  hiW   r40      h1   3.5
339 medP  hiW   r40      h1   3.5
340 medP  hiW   r40      h1   3.5
341 medP  hiW   r40      h1   3.5
342 medP  hiW   r40      h1   3.5
343 medP  hiW   r40      h2   5.0
344 medP  hiW   r40      h2   4.0
345 medP  hiW   r40      h2   2.0
346 medP  hiW   r40      h2   6.0
347 medP  hiW   r40      h2   1.0
348 medP  hiW   r40      h2   3.0
349 medP  hiW   r50      h1   3.5
350 medP  hiW   r50      h1   3.5
351 medP  hiW   r50      h1   3.5
352 medP  hiW   r50      h1   3.5
353 medP  hiW   r50      h1   3.5
354 medP  hiW   r50      h1   3.5
355 medP  hiW   r50      h2   6.0
356 medP  hiW   r50      h2   2.0
357 medP  hiW   r50      h2   5.0
358 medP  hiW   r50      h2   1.0
359 medP  hiW   r50      h2   3.0
360 medP  hiW   r50      h2   4.0
361 medP  loW   r10      h1   5.0
362 medP  loW   r10      h1   4.0
363 medP  loW   r10      h1   6.0
364 medP  loW   r10      h1   1.0
365 medP  loW   r10      h1   3.0
366 medP  loW   r10      h1   2.0
367 medP  loW   r10      h2   5.0
368 medP  loW   r10      h2   6.0
369 medP  loW   r10      h2   1.0
370 medP  loW   r10      h2   4.0
371 medP  loW   r10      h2   3.0
372 medP  loW   r10      h2   2.0
373 medP  loW   r20      h1   4.0
374 medP  loW   r20      h1   6.0
375 medP  loW   r20      h1   1.0
376 medP  loW   r20      h1   3.0
377 medP  loW   r20      h1   5.0
378 medP  loW   r20      h1   2.0
379 medP  loW   r20      h2   4.0
380 medP  loW   r20      h2   2.0
381 medP  loW   r20      h2   1.0
382 medP  loW   r20      h2   6.0
383 medP  loW   r20      h2   3.0
384 medP  loW   r20      h2   5.0
385 medP  loW   r30      h1   1.0
386 medP  loW   r30      h1   3.0
387 medP  loW   r30      h1   5.0
388 medP  loW   r30      h1   2.0
389 medP  loW   r30      h1   6.0
390 medP  loW   r30      h1   4.0
391 medP  loW   r30      h2   5.0
392 medP  loW   r30      h2   4.0
393 medP  loW   r30      h2   2.0
394 medP  loW   r30      h2   1.0
395 medP  loW   r30      h2   6.0
396 medP  loW   r30      h2   3.0
397 medP  loW   r40      h1   3.5
398 medP  loW   r40      h1   3.5
399 medP  loW   r40      h1   3.5
400 medP  loW   r40      h1   3.5
401 medP  loW   r40      h1   3.5
402 medP  loW   r40      h1   3.5
403 medP  loW   r40      h2   1.0
404 medP  loW   r40      h2   2.5
405 medP  loW   r40      h2   4.0
406 medP  loW   r40      h2   6.0
407 medP  loW   r40      h2   2.5
408 medP  loW   r40      h2   5.0
409 medP  loW   r50      h1   3.5
410 medP  loW   r50      h1   3.5
411 medP  loW   r50      h1   3.5
412 medP  loW   r50      h1   3.5
413 medP  loW   r50      h1   3.5
414 medP  loW   r50      h1   3.5
415 medP  loW   r50      h2   2.0
416 medP  loW   r50      h2   6.0
417 medP  loW   r50      h2   3.0
418 medP  loW   r50      h2   4.0
419 medP  loW   r50      h2   5.0
420 medP  loW   r50      h2   1.0
421 medP MedW   r10      h1   1.0
422 medP MedW   r10      h1   2.0
423 medP MedW   r10      h1   3.0
424 medP MedW   r10      h2   3.0
425 medP MedW   r10      h2   2.0
426 medP MedW   r10      h2   1.0
427 medP MedW   r20      h1   2.0
428 medP MedW   r20      h1   1.0
429 medP MedW   r20      h1   3.0
430 medP MedW   r20      h2   2.0
431 medP MedW   r20      h2   1.0
432 medP MedW   r20      h2   3.0
433 medP MedW   r30      h1   2.0
434 medP MedW   r30      h1   1.0
435 medP MedW   r30      h1   3.0
436 medP MedW   r30      h2   1.0
437 medP MedW   r30      h2   3.0
438 medP MedW   r30      h2   2.0
439 medP MedW   r40      h1   2.0
440 medP MedW   r40      h1   2.0
441 medP MedW   r40      h1   2.0
442 medP MedW   r40      h2   3.0
443 medP MedW   r40      h2   1.0
444 medP MedW   r40      h2   2.0
445 medP MedW   r50      h1   2.0
446 medP MedW   r50      h1   2.0
447 medP MedW   r50      h1   2.0
448 medP MedW   r50      h2   1.0
449 medP MedW   r50      h2   2.0
450 medP MedW   r50      h2   3.0


# exercice
# Calculate the average proportion of roots in each layer for each P/W treatment combination


# fitting models to different groups
# first, sum all root layers and harvests per pot
totalroot <- ddply(rootdata.melt, .(P, W, ID), numcolwise(sum))

# function to perform linear model
mymod <- function(x) {
  mod <- lm(value ~ P, data = x)
  anova(mod)
}
# use ddply
allmods <- dlply(totalroot, .(W), mymod)