The R Project SVN R-packages

Rev

Rev 5789 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 5789 Rev 5793
Line 105... Line 105...
105 
plot(air.boot)
 105 
plot(air.boot)
106 
 
 106 
 
107 
# In the difference of means example for the last two series of the
107 
# In the difference of means example for the last two series of the
108 
# gravity data
 108 
# gravity data
109 
grav1 <- gravity[as.numeric(gravity[, 2]) >= 7, ]
 109 
grav1 <- gravity[as.numeric(gravity[, 2]) >= 7, ]
110 
grav.fun <- function(dat, w)
 110 
grav.fun <- function(dat, w) {
111 
{    strata <- tapply(dat[, 2], as.numeric(dat[, 2]))
 111 
     strata <- tapply(dat[, 2], as.numeric(dat[, 2]))
112 
     d <- dat[, 1]
 112 
     d <- dat[, 1]
113 
     ns <- tabulate(strata)
 113 
     ns <- tabulate(strata)
114 
     w <- w/tapply(w, strata, sum)[strata]
 114 
     w <- w/tapply(w, strata, sum)[strata]
115 
     mns <- tapply(d * w, strata, sum)
 115 
     mns <- as.vector(tapply(d * w, strata, sum)) # drop names
116 
     mn2 <- tapply(d * d * w, strata, sum)
 116 
     mn2 <- tapply(d * d * w, strata, sum)
117 
     s2hat <- sum((mn2 - mns^2)/ns)
 117 
     s2hat <- sum((mn2 - mns^2)/ns)
118 
     c(mns[2]-mns[1], s2hat)
 118 
     c(mns[2] - mns[1], s2hat)
119 
}
 119 
}
120 
 
 120 
 
121 
grav.boot <- boot(grav1, grav.fun, R = 499, stype = "w", strata = grav1[, 2])
 121 
grav.boot <- boot(grav1, grav.fun, R = 499, stype = "w", strata = grav1[, 2])
122 
plot(grav.boot)
 122 
plot(grav.boot)
123 
# now suppose we want to look at the studentized differences.
 123 
# now suppose we want to look at the studentized differences.
124 
grav.z <- (grav.boot$t[, 1]-grav.boot$t0[1])/sqrt(grav.boot$t[, 2])
 124 
grav.z <- (grav.boot$t[, 1]-grav.boot$t0[1])/sqrt(grav.boot$t[, 2])
125 
plot(grav.boot, t = grav.z, t0 = 0)
 125 
plot(grav.boot, t = grav.z, t0 = 0)
126 
 
 126 
 
127 
# In this example we look at the one of the partial correlations for the
 127 
# In this example we look at the one of the partial correlations for the
128 
# head dimensions in the dataset frets.
 128 
# head dimensions in the dataset frets.
129 
pcorr <- function( x )
 129 
pcorr <- function(x) {
130 
{
- 
 
131 
#  Function to find the correlations and partial correlations between
 130 
#  Function to find the correlations and partial correlations between
132 
#  the four measurements.
 131 
#  the four measurements.
133 
     v <- cor(x);
 132 
     v <- cor(x)
134 
     v.d <- diag(var(x));
 133 
     v.d <- diag(var(x))
135 
     iv <- solve(v);
 134 
     iv <- solve(v)
136 
     iv.d <- sqrt(diag(iv));
 135 
     iv.d <- sqrt(diag(iv))
137 
     iv <- - diag(1/iv.d) \%*\% iv \%*\% diag(1/iv.d);
 136 
     iv <- - diag(1/iv.d) \%*\% iv \%*\% diag(1/iv.d)
138 
     q <- NULL;
137 
     q <- NULL
139 
     n <- nrow(v);
 138 
     n <- nrow(v)
140 
     for (i in 1:(n-1))
139 
     for (i in 1:(n-1))
141 
          q <- rbind( q, c(v[i, 1:i], iv[i, (i+1):n]) );
 140 
          q <- rbind( q, c(v[i, 1:i], iv[i, (i+1):n]) )
142 
     q <- rbind( q, v[n,] );
 141 
     q <- rbind( q, v[n, ] )
143 
     diag(q) <- round(diag(q));
 142 
     diag(q) <- round(diag(q))
144 
     q
 143 
     q
145 
}
 144 
}
146 
 
 145 
 
147 
frets.fun <- function(data, i) {
 146 
frets.fun <- function(data, i) {
148 
     d <- data[i,];
 147 
     d <- data[i, ]
149 
     v <- pcorr( d );
 148 
     v <- pcorr(d)
150 
     c(v[1, ], v[2, ], v[3, ], v[4, ])
 149 
     c(v[1, ], v[2, ], v[3, ], v[4, ])
151 
}
 150 
}
152 
frets.boot <- boot(log(as.matrix(frets)),  frets.fun,  R = 999)
 151 
frets.boot <- boot(log(as.matrix(frets)),  frets.fun,  R = 999)
153 
plot(frets.boot, index = 7, jack = TRUE, stinf = FALSE, useJ = FALSE)
 152 
plot(frets.boot, index = 7, jack = TRUE, stinf = FALSE, useJ = FALSE)
154 
}
 153 
}