has at least a reasonably small P-value of 0.088. This
is because the significant difference does not involve the baseline:
it’s actually between flat
with a positive slope and
upsloping
with a negative one.
\(\blacksquare\)
- We are going to make a large number of predictions. Create and save
a data frame that contains predictions for all combinations of representative
values for all the variables in the model that came out of
step
. By “representative” I mean all the values for a categorical variable, and the five-number summary for a numeric variable. (Note that you will get a lot of predictions.)
Solution
This is exactly what predictions
with variables
does, so that you don’t have to do anything to get hold of the representative variables yourself.
The hard work is in listing all the variables. The easiest way to make sure you have them all is to look at the summary
of your best model (mine was called heart.3
) first, and copy them from the Call
at the top. This is easier than looking at the table of Coefficients (or tidy
output) because for categorical variables like pain.type
you will have to distinguish the name of the variable from its levels. For example, the table of Coefficients has pain.typeatypical
and pain.typenonanginal
.  That is fine — as long as you don’t display them all for a grader to have to page through!
\(\blacksquare\)
- Find the largest predicted probability (which is the predicted probability of heart disease) and display all the variables that it was a prediction for.
Solution
The (at current writing) approved way to do this is to use slice_max
. This finds the rows with maximum value(s) on a variable, which is exactly what we want.
Variations:
- by using
prop
instead ofn
, you can display theproportion
of rows with the highest values on your variable, such as the 10% of rows with the highest predicted probabilities withprop = 0.10
- there is also
slice_min
that displays the rows with the lowest values on a variable, or the input proportion of rows with the lowest values - there are alse
slice_head
andslice_tail
that display the first and last (respectively) rows in a dataframe. The default display of a dataframe in an R Notebook is thus
%>% slice_head(n = 10) p
## rowid type predicted std.error conf.low conf.high heart.disease sex
## 1 1 response 0.4840884 0.26015315 0.10858288 0.8784638 no male
## 2 2 response 0.7880319 0.17215432 0.33028444 0.9655472 no male
## 3 3 response 0.4043889 0.26993345 0.07019033 0.8592830 no male
## 4 4 response 0.7447668 0.19534420 0.28024364 0.9562714 no male
## 5 5 response 0.9203905 0.07531898 0.60657581 0.9885967 no male
## 6 6 response 0.6786011 0.23916970 0.19750828 0.9476802 no male
## 7 7 response 0.9657764 0.03922606 0.73379447 0.9965505 no male
## 8 8 response 0.9911354 0.01060165 0.91307661 0.9991604 no male
## 9 9 response 0.9533128 0.05461959 0.64821589 0.9956000 no male
## 10 10 response 0.2517945 0.20668885 0.03771116 0.7429256 no male
## pain.type resting.bp serum.chol max.hr oldpeak slope colored
## 1 asymptomatic 94 126 71 0 flat 0
## 2 asymptomatic 94 126 71 0 flat 0
## 3 asymptomatic 94 126 71 0 flat 0
## 4 asymptomatic 94 126 71 0 flat 1
## 5 asymptomatic 94 126 71 0 flat 1
## 6 asymptomatic 94 126 71 0 flat 1
## 7 asymptomatic 94 126 71 0 flat 3
## 8 asymptomatic 94 126 71 0 flat 3
## 9 asymptomatic 94 126 71 0 flat 3
## 10 asymptomatic 94 126 71 0 upsloping 0
## thal
## 1 normal
## 2 reversible
## 3 fixed
## 4 normal
## 5 reversible
## 6 fixed
## 7 normal
## 8 reversible
## 9 fixed
## 10 normal
except that the default display also tells you how many rows there are altogether.
- you may have run into
slice_sample
, which displays a randomly-chosen number or proportion of rows from a dataframe. This is useful after you read in a dataframe from a file, if you want to get a sense of what kind of values you have in your dataframe (for example, if they are ordered by something and looking at the first ten rows won’t tell you the whole story, such as having males listed first and you want to check that there are some females as well):
%>% slice_sample(n = 10) p
## rowid type predicted std.error conf.low conf.high heart.disease
## 1 53751 response 0.7256533 0.376720147 0.06087308 0.9908201 no
## 2 47026 response 0.2036584 0.238000424 0.01420553 0.8194534 no
## 3 19771 response 0.6917220 0.197561555 0.26743943 0.9323914 no
## 4 45758 response 0.6197210 0.393948586 0.05797736 0.9773505 no
## 5 26121 response 0.8578826 0.205782593 0.18089296 0.9939759 no
## 6 76988 response 0.3835134 0.357924015 0.03101496 0.9236112 no
## 7 5069 response 0.9971035 0.005399105 0.89819727 0.9999256 no
## 8 5995 response 0.8978355 0.086041059 0.58295851 0.9822224 no
## 9 33785 response 0.9983106 0.002701200 0.96240473 0.9999267 no
## 10 86160 response 0.8121677 0.297026385 0.08690634 0.9949350 no
## sex pain.type resting.bp serum.chol max.hr oldpeak slope colored
## 1 male typical 200 564 153.5 1.6 downsloping 0
## 2 male typical 130 245 71.0 0.8 downsloping 0
## 3 male nonanginal 130 213 166.0 0.0 flat 3
## 4 male typical 120 564 166.0 1.6 downsloping 0
## 5 male nonanginal 200 281 133.0 6.2 upsloping 0
## 6 female nonanginal 140 245 153.5 6.2 upsloping 0
## 7 male asymptomatic 120 564 133.0 6.2 downsloping 0
## 8 male asymptomatic 130 213 71.0 1.6 flat 0
## 9 male atypical 130 245 153.5 6.2 flat 3
## 10 female atypical 120 564 153.5 6.2 flat 0
## thal
## 1 fixed
## 2 normal
## 3 normal
## 4 reversible
## 5 fixed
## 6 reversible
## 7 reversible
## 8 normal
## 9 reversible
## 10 fixed
If you didn’t think of slice_max
, there are lots of other ways. Find one.  Are they consistent?
Solution
Since we were predicting the probability of heart disease, a
more positive slope in the model from step
will be
associated with a higher probability of heart disease. So,
there, we are looking for a couple of things: if the variable
is a factor, we’re looking for the level with the most
positive slope (bearing in mind that this might be the
baseline), and for a numeric variable, if the slope is
positive, a high value is associated with heart
disease, and if negative, a low value.
Bearing that in mind, we go back to my
summary(heart.3)
:
summary(heart.3)
##
## Call:
## glm(formula = factor(heart.disease) ~ sex + pain.type + resting.bp +
## serum.chol + max.hr + oldpeak + slope + colored + thal, family = "binomial",
## data = heart)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7103 -0.4546 -0.1442 0.3864 2.7121
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.818418 2.550437 -1.889 0.058858 .
## sexmale 1.850559 0.561583 3.295 0.000983 ***
## pain.typeatypical -1.268233 0.604488 -2.098 0.035903 *
## pain.typenonanginal -2.086204 0.486591 -4.287 1.81e-05 ***
## pain.typetypical -2.532340 0.748941 -3.381 0.000722 ***
## resting.bp 0.024125 0.011077 2.178 0.029410 *
## serum.chol 0.007142 0.003941 1.812 0.069966 .
## max.hr -0.020373 0.010585 -1.925 0.054262 .
## oldpeak 0.467028 0.233280 2.002 0.045284 *
## slopeflat 0.859564 0.922749 0.932 0.351582
## slopeupsloping -0.165832 0.991474 -0.167 0.867167
## colored 1.134561 0.261547 4.338 1.44e-05 ***
## thalnormal 0.323543 0.813442 0.398 0.690818
## thalreversible 1.700314 0.805127 2.112 0.034699 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 370.96 on 269 degrees of freedom
## Residual deviance: 173.57 on 256 degrees of freedom
## AIC: 201.57
##
## Number of Fisher Scoring iterations: 6
sex
: being male has the higher risk, by a lotpain
: all the slopes shown are negative, so the highest risk goes with the baseline oneasymptomatic
.resting.bp
: positive slope, so higher risk with higher value.serum.chol
: same.max.hr
: negative slope, so greatest risk with smaller value.oldpeak
: positive slope, greater risk with higher value again.slope
:flat
has greatest risk.colored
: positive slope, so beware of higher value.thal
:reversible
has greatest risk.
Then we can do the same thing for the prediction.
and the highest prediction:
%>% slice_max(predicted, n = 1) p
## rowid type predicted std.error conf.low conf.high heart.disease sex
## 1 13049 response 0.9999983 3.918393e-06 0.999837 1 no male
## pain.type resting.bp serum.chol max.hr oldpeak slope colored thal
## 1 asymptomatic 200 564 71 6.2 flat 3 reversible
For the numerical variables, we may need to check back to the previous part to see whether the value shown was high or low. Once you have done that, you can see that the variable values for the highest predicted probability do indeed match the ones we thought should be the highest risk.
Extra: the interesting thing about this is that even after adjusting for all of the other variables, there is a greater risk of heart disease if you are male (and the model shows that the risk is much greater). That is to say, it’s being male that makes the difference, not the fact that any of the other variables are different for males.
It’s rather difficult to scan 108,000 predictions to see the effect of being male, but we can do this:
predictions(heart.3, variables = "sex")
## rowid type predicted std.error heart.disease pain.type resting.bp
## 1 1 response 0.5849307 0.1111510 no asymptomatic 131.3444
## 2 2 response 0.1813080 0.0810187 no asymptomatic 131.3444
## serum.chol max.hr oldpeak slope colored thal sex
## 1 249.6593 149.6778 1.05 upsloping 0.6703704 normal male
## 2 249.6593 149.6778 1.05 upsloping 0.6703704 normal female
What this does is to choose a single representative value for all the other variables: the mean for a quantitative variable like resting blood pressure, and the most common category for a categorical variable like pain.type
. If you scan all the way along the two rows, you find that the values for all the variables are the same except for sex
at the end. The predicted probabilities of heart disease are very different for males and females (much higher for males), especially given that all else really is equal.
To see this graphically, we can use plot_cap
, and we can include another variable such as resting blood pressure. It’s better to list the quantitative one first:
plot_cap(heart.3, condition = c("resting.bp", "sex"))
As the resting blood pressure increases, the probability of heart disease increases, but the blue line for males is well above the red one for females all the way across. For example, for a 50% chance of heart disease, this will happen for males with a resting blood pressure of about 120, but for females not until the resting blood pressure reaches 190!
Perhaps, therefore, the easiest way to avoid a heart attack is to not be male!
\(\blacksquare\)
26.17 Successful breastfeeding
A regular pregnancy lasts 40 weeks, and a baby that is born at or before 33 weeks is called “premature”. The number of weeks at which a baby is born is called its “gestational age”. Premature babies are usually smaller than normal and may require special care. It is a good sign if, when the mother and baby leave the hospital to go home, the baby is successfully breastfeeding.
The data in link are from a study of 64 premature infants. There are three columns: the gestational age (a whole number of weeks), the number of babies of that gestational age that were successfully breastfeeding when they left the hospital, and the number that were not. (There were multiple babies of the same gestational age, so the 64 babies are summarized in 6 rows of data.)
- Read the data into R and display the data frame.
Solution
No great challenge here, I hope:
<- "http://ritsokiguess.site/datafiles/breastfeed.csv"
my_url <- read_csv(my_url) breastfeed
## Rows: 6 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (3): gest.age, bf.yes, bf.no
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
breastfeed
## # A tibble: 6 × 3
## gest.age bf.yes bf.no
## <dbl> <dbl> <dbl>
## 1 28 2 4
## 2 29 2 3
## 3 30 7 2
## 4 31 7 2
## 5 32 16 4
## 6 33 14 1
That looks reasonable.
\(\blacksquare\)
- Verify that there were indeed 64 infants, by having R do a suitable calculation on your data frame that gives the right answer for the right reason.
Solution
The second and third columns are all frequencies, so it’s a question of adding them up. For example:
%>% summarize(total = sum(bf.yes) + sum(bf.no)) breastfeed
## # A tibble: 1 × 1
## total
## <dbl>
## 1 64
or if you want to go nuts (this one pivot-longers all the frequencies together into one column and then adds them up):
%>%
breastfeed pivot_longer(bf.yes:bf.no, names_to="yesno", values_to="freq") %>%
summarize(total = sum(freq))
## # A tibble: 1 × 1
## total
## <dbl>
## 1 64
Find a way to get it done. If it works and it does the right thing, it’s good.
Do not copy the numbers out of your data frame, type them in again and use R to add them up. Do something with your data frame as you read it in.
\(\blacksquare\)
- Do you think, looking at the data, that there is a relationship between gestational age and whether or not the baby was successfully breastfeeding when it left the hospital? Explain briefly.
Solution
The babies with the youngest gestational age (the most premature) were mostly not breastfeeding when they left the hospital. Most of the 30- and 31-week babies were breastfeeding, and almost all of the 32- and 33-week babies were breastfeeding. So I think there will be a relationship: as gestational age increases, the probability that the baby will be breastfeeding will also increase. (This, looking ahead, suggests a positive slope in a logistic regression.)
\(\blacksquare\)
- Why is logistic regression a sensible technique to use here? Explain briefly.
Solution
The response variable is a yes/no: whether or not an infant is breastfeeding. We want to predict the probability of the response being in one or the other category. This is what logistic regression does. (The explanatory variable(s) are usually numerical, as here, but they could be factors as well, or a mixture. The key is the kind of response. The number of babies that are successfully breastfeeding at a certain gestational age is modelled as binomial with \(n\) being the total number of babies of that gestational age, and \(p\) being something that might depend, and here does depend, on gestational age.)
\(\blacksquare\)
- Fit a logistic regression to predict the probability that an infant will be breastfeeding from its gestational age. Show the output from your logistic regression.
Solution
These are summarized data, rather than one infant per line, so
what we have to do is to make a two-column response “matrix”,
successes in the first column and failures in the second, and then
predict that from gestational age. (That’s why this was
three marks rather than two.)
So, let’s make the response
first:
<- with(breastfeed, cbind(bf.yes, bf.no))
response response
## bf.yes bf.no
## [1,] 2 4
## [2,] 2 3
## [3,] 7 2
## [4,] 7 2
## [5,] 16 4
## [6,] 14 1
or, more Tidyverse-like, but we have to remember to turn it into a
matrix
:
<- breastfeed %>%
response select(starts_with("bf")) %>%
as.matrix()
response
## bf.yes bf.no
## [1,] 2 4
## [2,] 2 3
## [3,] 7 2
## [4,] 7 2
## [5,] 16 4
## [6,] 14 1
I used a select-helper, because what immediately came to me was that
the names of the columns I wanted started with bf
, but
whatever way you have that works is good.
Now we fit the logistic regression:
.1 <- glm(response ~ gest.age, data = breastfeed, family = "binomial")
breastfeedsummary(breastfeed.1)
##
## Call:
## glm(formula = response ~ gest.age, family = "binomial", data = breastfeed)
##
## Deviance Residuals:
## 1 2 3 4 5 6
## -0.1472 -0.4602 0.8779 0.1114 -0.6119 0.3251
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -16.7198 6.0630 -2.758 0.00582 **
## gest.age 0.5769 0.1977 2.918 0.00352 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 11.1772 on 5 degrees of freedom
## Residual deviance: 1.4968 on 4 degrees of freedom
## AIC: 19.556
##
## Number of Fisher Scoring iterations: 4
\(\blacksquare\)
- Does the significance or non-significance of the slope of
gest.age
surprise you? Explain briefly.
Solution
The slope is significant (P-value 0.0035 is much less than 0.05). We said above that we expected there to be a relationship between gestational age and whether or not the baby was breastfeeding, and this significant slope is confirming that there is a relationship. So this is exactly what we expected to see, and not a surprise at all. If you concluded above that you did not see a relationship, you should colour yourself surprised here. Consistency.
\(\blacksquare\)
- Is your slope (in the
Estimate
column) positive or negative? What does that mean, in terms of gestational ages and breastfeeding? Explain briefly.
Solution
My slope is 0.5769, positive. That means that as the explanatory variable, gestational age, increases, the probability of the event (that the baby is breastfeeding) also increases. This is also what I observed above: almost all of the near-term (large gestational age) babies were breastfeeding, whereas a fair few of the small-gestational-age (very premature) ones were not.
We know that the event is “is breastfeeding” rather than “is not”, because the first column of our response matrix is “yes” rather than “no”:
response
## bf.yes bf.no
## [1,] 2 4
## [2,] 2 3
## [3,] 7 2
## [4,] 7 2
## [5,] 16 4
## [6,] 14 1
(If you had happened to make your response matrix the other way around, the event would have been “is not breastfeeding”, and then your slope would have been the same size but negative.)
\(\blacksquare\)
- Obtain the predicted probabilities that an infant will successfully breastfeed for a representative collection of gestational ages.
Solution
This is the easier version of predictions
, where you do not
have to create a data frame of values to predict from.
Thus, you only need something like this, with variables
(and no datagrid
first):
predictions(breastfeed.1, variables = "gest.age")
## rowid type predicted std.error conf.low conf.high gest.age
## 1 1 response 0.3620282 0.14035792 0.1470846 0.6512432 28.0
## 2 2 response 0.5025820 0.11259364 0.2947524 0.7095232 29.0
## 3 3 response 0.7059259 0.06554832 0.5638536 0.8167588 30.5
## 4 4 response 0.8508177 0.05270073 0.7165154 0.9278964 32.0
## 5 5 response 0.9103511 0.04606209 0.7706029 0.9684505 33.0
You can see that the predicted probabilities go steadily up as the gestational age goes up, just as we would have expected.
Extra: if you only wanted certain gestational ages, for example 25, 30 and 35, you would do that like this:
<- datagrid(model = breastfeed.1, gest.age = c(25, 30, 35))
new predictions(breastfeed.1, newdata = new)
## rowid type predicted std.error conf.low conf.high gest.age
## 1 1 response 0.09134905 0.09599845 0.01031219 0.4923797 25
## 2 2 response 0.64272897 0.07771617 0.48098221 0.7773976 30
## 3 3 response 0.96987261 0.02686505 0.84153715 0.9949018 35
Or, if you wanted to make a graph of the observed and predicted proportions/probabilities, you could do something like this:
plot_cap(breastfeed.1, condition = "gest.age")
or, if you wanted to add the data to this:
%>%
breastfeed mutate(total = bf.yes + bf.no,
proportion = bf.yes / total) -> d
plot_cap(breastfeed.1, condition = "gest.age") +
geom_point(data = d, aes(x = gest.age, y = proportion, size = total), inherit.aes = FALSE)
What did I do there? I first created some new variables:
total
is the total number of babies of each gestational age, and
proportion
is the observed proportion of breastfeeding babies at
each gestational age (number of yes divided by total).
pred
which are the predictions we did above. Then I repeat the plot_cap
, and on that plot
I add the observed proportions against gestational age
(as points). To do that, I need to change which dataframe I am using (the temporary one called d
with the proportions and totals in it), and the variables I am plotting (which are now the ones in d
).
The way I
do that is to put an aes
inside the geom_line
to
say “use this x
and y
instead”. I also
wanted to draw attention to the gestational ages where more babies
were observed; I did this by making the size of the plotted
points proportional to how many babies there were at that gestational
age (which was the quantity total
I calculated above).
The legend for total
tells you what size point
corresponds to how many total babies. The final thing is that plot_cap
defines some other things (like the min and max of the confidence interval) that we would also have to supply values for were it not for the inherit.aes
at the end; this means “don’t use anything from the original ggplot
(the plot_cap
), but supply everything yourself”.
The idea is that the observed and predicted should be reasonably close, or at least not grossly different, especially when there is a lot of data (the circles are large), and I think they are close, which indicates that our model is doing a good job. You can see that there is not much data on the left, with small gestational ages, so the confidence interval around the predictions is wider there. On the right, where there is lots of data, the interval is narrower. When the gestational age is large (the baby is closer to being full term rather than premature), there is a good chance that the baby will be able to breastfeed, and we are fairly sure about that.
\(\blacksquare\)
26.18 Making it over the mountains
In 1846, the Donner party (Donner and Reed families) left Springfield, Illinois for California in covered wagons. After reaching Fort Bridger, Wyoming, the leaders decided to find a new route to Sacramento. They became stranded in the eastern Sierra Nevada mountains at a place now called Donner Pass, when the region was hit by heavy snows in late October. By the time the survivors were rescued on April 21, 1847, 40 out of 87 had died.
After the rescue, the age and gender of each person in the party was recorded, along with whether they survived or not. The data are in link.
- Read in the data and display its structure. Does that agree with the description in the previous paragraph?
Solution
Nothing very new here:
<- "http://ritsokiguess.site/datafiles/donner.txt"
my_url <- read_delim(my_url, " ") donner
## Rows: 45 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: " "
## chr (2): gender, survived
## dbl (1): age
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
donner
## # A tibble: 45 × 3
## age gender survived
## <dbl> <chr> <chr>
## 1 23 male no
## 2 40 female yes
## 3 40 male yes
## 4 30 male no
## 5 28 male no
## 6 40 male no
## 7 45 female no
## 8 62 male no
## 9 65 male no
## 10 45 female no
## # … with 35 more rows
The ages look like ages, and the other variables are categorical as promised, with the right levels. So this looks sensible.
I usually write (or rewrite) the data description myself, but the one I found for these data at link was so nice that I used it as is. You can see that the data format used on that website is not as nice as ours. (I did some work to make it look nicer for you, with proper named categories rather than 0 and 1.)
\(\blacksquare\)
- Make graphical or numerical summaries for each pair of
variables. That is, you should make a graph or numerical summary for
each of
age
vs.gender
,age
vs.
survived
andgender
vs.survived
. In choosing the kind of graph or summary that you will use, bear in mind thatsurvived
andgender
are factors with two levels.
Solution
Thinking about graphs first: we have two numeric-vs-factor graphs
(the ones involving age
), which I think should be
boxplots, though they could be something like side-by-side
histograms (if you are willing to grapple with
facet_grid
). The other two variables are both
factors, so a good graph for them would be something
like a grouped bar plot (as in the question about parasites earlier).
If you prefer numerical summaries: you could do mean age (or some
other summary of age) for each group defined by gender or
survivorship, and you could do a cross-tabulation of frequencies
for gender and survival. I’ll take any mixture of graphs and
numerical summaries that addresses each pair of variables somehow
and summarizes them in a sensible way.
Starting with age
vs. gender
:
ggplot(donner, aes(x = gender, y = age)) + geom_boxplot()
or:
ggplot(donner, aes(x = age)) + geom_histogram(bins = 10) + facet_grid(gender ~ .)
or:
%>% group_by(gender) %>% summarize(m = mean(age)) donner
## # A tibble: 2 × 2
## gender m
## <chr> <dbl>
## 1 female 31.1
## 2 male 32.2
Age vs. survived
is the same idea:
ggplot(donner, aes(x = survived, y = age)) + geom_boxplot()
or:
ggplot(donner, aes(x = age)) + geom_histogram(bins = 10) + facet_grid(survived ~ .)
or:
%>% group_by(survived) %>% summarize(m = mean(age)) donner
## # A tibble: 2 × 2
## survived m
## <chr> <dbl>
## 1 no 35.5
## 2 yes 27.2
For survived
against gender
, the obvious thing is a
cross-tabulation, gotten like this:
with(donner, table(gender, survived))
## survived
## gender no yes
## female 5 10
## male 20 10
or like this:
%>% group_by(gender, survived) %>% summarize(n = n()) donner
## `summarise()` has grouped output by 'gender'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 3
## # Groups: gender [2]
## gender survived n
## <chr> <chr> <int>
## 1 female no 5
## 2 female yes 10
## 3 male no 20
## 4 male yes 10
For a graph, borrow the grouped bar-plot idea from the parasites question:
ggplot(donner, aes(x = gender, fill = survived)) + geom_bar(position = "dodge")
I think this way around for x
and fill
is better,
since we want to ask something like “out of the females, how many survived”
(that is, gender is explanatory and survival response).
\(\blacksquare\)
- For each of the three graphs or summaries in the previous question, what do they tell you about the relationship between the pair of variables concerned? Explain briefly.
Solution
I think the graphs are the easiest thing to interpret, but use whatever you got:
The females on average were younger than the males. (This was still true with the medians, even though those very old males might have pulled the mean up.)
The people who survived were on average younger than those who didn’t (or, the older people tended not to survive).
A greater proportion of females survived than males.
A relevant point about each relationship is good.
\(\blacksquare\)
- Fit a logistic regression predicting survival from age and gender. Display the summary.
Solution
Each row of the data frame is one person, so we can use the
survived
column without going through that two-column
response business.
However, the response variable survived
is a categorical
variable expressed as text, so we need to make it into a
factor
first. Either create a new variable that is the
factor version of survived
, or do it right in the glm
:
.1 <- glm(factor(survived) ~ age + gender, family = "binomial", data = donner)
donnersummary(donner.1)
##
## Call:
## glm(formula = factor(survived) ~ age + gender, family = "binomial",
## data = donner)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7445 -1.0441 -0.3029 0.8877 2.0472
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.23041 1.38686 2.329 0.0198 *
## age -0.07820 0.03728 -2.097 0.0359 *
## gendermale -1.59729 0.75547 -2.114 0.0345 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 61.827 on 44 degrees of freedom
## Residual deviance: 51.256 on 42 degrees of freedom
## AIC: 57.256
##
## Number of Fisher Scoring iterations: 4
We ought to take a moment to think about what is being predicted here:
levels(factor(donner$survived))
## [1] "no" "yes"
The baseline is the first of these, no
, and the thing that is
predicted is the probability of the second one, yes
(that is,
the probability of surviving).
\(\blacksquare\)
- Do both explanatory variables have an impact on survival? Does that seem to be consistent with your numerical or graphical summaries? Explain briefly.
Solution
Both explanatory variables have a P-value (just) less than 0.05, so they both have an impact on survival: taking either one of them out of the logistic regression would be a mistake. To see if this makes sense, go back to your plots or summaries, the ones involving survival. For age, the mean or median age of the survivors was less than for the people who died, by five year (median) or eight years (mean), so it makes sense that there would be an age effect. For gender, two-thirds of the women survived and two-thirds of the men died, so there ought to be a gender effect and is.
\(\blacksquare\)
- Are the men typically older, younger or about the same age as
the women? Considering this, explain carefully what the negative
gendermale
slope in your logistic regression means.
Solution
The men are typically older than the women.
The negative (and significant) gendermale
slope means that
the probability of a male surviving is less than that of a woman
of the same age. Or, though the males are typically
older, even after you allow for that, the males have worse
survival. (genderfemale
is the baseline.)
You need to get at the idea of “all else equal” when you’re
assessing regression slopes of any kind: regular regression,
logistic regression or survival analysis (coming up later). That’s
why I said “carefully” in the question. If I say “carefully” or
“precisely”, a complete answer is looking for a specific thing to
show that you understand the issue completely.
\(\blacksquare\)
- Obtain predicted probabilities of survival for each combination of some representative ages and of the two genders in this dataset.
Solution
This is the simplest version of predictions
, using variables
to pick the representative values for you:
predictions(donner.1, variables = c("age", "gender"))
## rowid type predicted std.error conf.low conf.high survived age
## 1 1 response 0.61302822 0.14835669 0.317415947 0.8436698 no 15
## 2 2 response 0.88668846 0.09181091 0.566195670 0.9791302 no 15
## 3 3 response 0.43935567 0.10638552 0.251566044 0.6462789 no 24
## 4 4 response 0.79470387 0.11330089 0.498099889 0.9378846 no 24
## 5 5 response 0.36433597 0.09489750 0.204289024 0.5613180 no 28
## 6 6 response 0.73898501 0.12165366 0.451291689 0.9069418 no 28
## 7 7 response 0.18316607 0.09209645 0.062884867 0.4283504 no 40
## 8 8 response 0.52554048 0.15989796 0.239648863 0.7956146 no 40
## 9 9 response 0.03076455 0.04325543 0.001845139 0.3527587 no 65
## 10 10 response 0.13553852 0.15389121 0.011807792 0.6729187 no 65
## gender
## 1 male
## 2 female
## 3 male
## 4 female
## 5 male
## 6 female
## 7 male
## 8 female
## 9 male
## 10 female
(There are five ages chosen, something like a five-number summary, and two genders here, so ten predictions.)
These, remember, are predicted probabilities of surviving.
\(\blacksquare\)
- Do your predictions support your conclusions from earlier about the effects of age and gender? Explain briefly.
Solution
We said before that the probability of survival was lower if the age was higher. This is confirmed here: for example, look at the odd-numbered rows 1, 3, 5, 7, 9, which are all females of increasing ages; the probability of survival decreases. (Or look at males, in the even-numbered rows; the effect is the same.)
To see the effect of gender, look at two predictions of different genders but the same age (eg. rows 1 and 2). The female is always predicted to have the higher survival probability. This is also what we saw before. The effect of gender is substantial, but not strongly significant, because we only have 45 observations, not so many when all we know about each person is whether they survived or not. I wanted you to think about these different ways to understand the model, and to understand that they all say the same thing, in different ways (and thus you can look at whichever of them is most convenient or comprehensible). For the logistic and survival models, I find looking at predictions to be the easiest way to understand what the model is saying.
\(\blacksquare\)
26.19 Who needs the most intensive care?
The “APACHE II” is a scale for assessing patients who arrive in the intensive care unit (ICU) of a hospital. These are seriously ill patients who may die despite the ICU’s best attempts. APACHE stands for “Acute Physiology And Chronic Health Evaluation”.49 The scale score is calculated from several physiological measurements such as body temperature, heart rate and the Glasgow coma scale, as well as the patient’s age. The final result is a score between 0 and 71, with a higher score indicating more severe health issues. Is it true that a patient with a higher APACHE II score has a higher probability of dying?
Data from one hospital are in link. The columns are: the APACHE II score, the total number of patients who had that score, and the number of patients with that score who died.
- Read in and display the data (however much of it displays). Why are you convinced that have the right thing?
Solution
Data values separated by one space, so:
<- "http://ritsokiguess.site/datafiles/apache.txt"
my_url <- read_delim(my_url, " ") icu
## Rows: 38 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: " "
## dbl (3): apache, patients, deaths
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
icu
## # A tibble: 38 × 3
## apache patients deaths
## <dbl> <dbl> <dbl>
## 1 0 1 0
## 2 2 1 0
## 3 3 4 1
## 4 4 11 0
## 5 5 9 3
## 6 6 14 3
## 7 7 12 4
## 8 8 22 5
## 9 9 33 3
## 10 10 19 6
## # … with 28 more rows
I had to stop and think about what to call the data frame, since one
of the columns is called apache
.
Anyway, I appear to have an apache
score between 0 and something, a
number of patients and a number of deaths (that is no bigger than the
number of patients). If you check the original data, the apache
scores go up to 41 and are all the values except for a few near the
end, so it makes perfect sense that there would be 38 rows.
Basically, any comment here is good, as long as you make one and it has something to do with the data.
apache
scores could be as high as 71, but I imagine a patient would
have to be very ill to get a score anywhere near that high.
\(\blacksquare\)
- Does each row of the data frame relate to one patient or sometimes to more than one? Explain briefly.
Solution
Sometimes to more than one. The number in the patients
column says how many patients that line refers to: that is to say
(for example) the row where apache
equals 6 represents
all the patients whose apache
score was 6, however many
of them there were (14 in this case).
I had to be careful with the wording because the first two rows of
the data frame actually do refer to only one patient each
(who survived in both cases), but the later rows do refer to more
than one patient.
\(\blacksquare\)
- Explain why this is the kind of situation where you need a
two-column response, and create this response variable, bearing in
mind that I will (later) want you to estimate the probability of
dying, given the
apache
score.
Solution
This needs a two-column response precisely because each row represents (or could represent) more than one observation. The two columns are the number of observations referring to the event of interest (dying), and the number of observations where that didn’t happen (survived). We don’t actually have the numbers of survivals, but we can calculate these by subtracting from the numbers of patients (since a patient must have either lived or died):
<- icu %>%
response mutate(survivals = patients - deaths) %>%
select(deaths, survivals) %>%
as.matrix()
response
## deaths survivals
## [1,] 0 1
## [2,] 0 1
## [3,] 1 3
## [4,] 0 11
## [5,] 3 6
## [6,] 3 11
## [7,] 4 8
## [8,] 5 17
## [9,] 3 30
## [10,] 6 13
## [11,] 5 26
## [12,] 5 12
## [13,] 13 19
## [14,] 7 18
## [15,] 7 11
## [16,] 8 16
## [17,] 8 19
## [18,] 13 6
## [19,] 7 8
## [20,] 6 7
## [21,] 9 8
## [22,] 12 2
## [23,] 7 6
## [24,] 8 3
## [25,] 8 4
## [26,] 2 4
## [27,] 5 2
## [28,] 1 2
## [29,] 4 3
## [30,] 4 1
## [31,] 3 0
## [32,] 3 0
## [33,] 1 0
## [34,] 1 0
## [35,] 1 0
## [36,] 1 0
## [37,] 1 0
## [38,] 0 1
noting that the deaths column has to come first since that’s
what we want the probability of. It has to be a matrix
, so
as.matrix
is the final step. You can quickly check that the
two numbers in each row add up to the number of patients
for
that row.
Or do everything outside of the data frame:
<- with(icu, patients - deaths)
survivals <- with(icu, cbind(deaths, survivals))
resp resp
## deaths survivals
## [1,] 0 1
## [2,] 0 1
## [3,] 1 3
## [4,] 0 11
## [5,] 3 6
## [6,] 3 11
## [7,] 4 8
## [8,] 5 17
## [9,] 3 30
## [10,] 6 13
## [11,] 5 26
## [12,] 5 12
## [13,] 13 19
## [14,] 7 18
## [15,] 7 11
## [16,] 8 16
## [17,] 8 19
## [18,] 13 6
## [19,] 7 8
## [20,] 6 7
## [21,] 9 8
## [22,] 12 2
## [23,] 7 6
## [24,] 8 3
## [25,] 8 4
## [26,] 2 4
## [27,] 5 2
## [28,] 1 2
## [29,] 4 3
## [30,] 4 1
## [31,] 3 0
## [32,] 3 0
## [33,] 1 0
## [34,] 1 0
## [35,] 1 0
## [36,] 1 0
## [37,] 1 0
## [38,] 0 1
class(resp)
## [1] "matrix" "array"
Or use the dollar sign instead of the with
s. Any of those is
good.
I have no objection to your displaying the response matrix.
\(\blacksquare\)
- Fit a logistic regression to estimate the probability of
death from the
apache
score, and display the results.
Solution
.1 <- glm(response ~ apache, family = "binomial", data = icu)
apachesummary(apache.1)
##
## Call:
## glm(formula = response ~ apache, family = "binomial", data = icu)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2508 -0.5662 0.1710 0.6649 2.3695
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.2903 0.2765 -8.282 < 2e-16 ***
## apache 0.1156 0.0160 7.227 4.94e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 106.009 on 37 degrees of freedom
## Residual deviance: 43.999 on 36 degrees of freedom
## AIC: 125.87
##
## Number of Fisher Scoring iterations: 4
My naming convention has gotten messed up again. This should really be
called deaths.1
or something like that, but that would be a
really depressing name.
\(\blacksquare\)
- Is there a significant effect of
apache
score on the probability of survival? Explain briefly.
Solution
A gimme two points. The P-value for apache
is \(4.94 \times 10^{-13}\), very small, so apache
score definitely
has an effect on the probability of survival.
\(\blacksquare\)
- Is the effect of a larger
apache
score to increase or to decrease the probability of death? Explain briefly.
Solution
The slope coefficient for apache
is 0.1156, positive, and
since we are modelling the probability of death (the first column
of the response matrix), this says that as apache
goes
up, the probability of death goes up as well.
If you made your response matrix with the columns the other way
around, the slope coefficient for apache
should be
\(-0.1156\), but the explanation should come to the same place,
because this says that the probability of survival goes down as
apache
goes up.
\(\blacksquare\)
- Obtain the predicted probability of death for a representative collection of the
apache
scores that were in the data set. Do your predictions behave as you would expect (from your earlier work)?
Solution
“Representative” is a clue that you can let predictions
pick the exact values using variables
:
predictions(apache.1, variables = "apache")
## rowid type predicted std.error conf.low conf.high apache
## 1 1 response 0.09192721 0.02308364 0.05560281 0.1482565 0.0
## 2 2 response 0.24341350 0.02605492 0.19602193 0.2980157 10.0
## 3 3 response 0.49110180 0.02936804 0.43391128 0.5485262 19.5
## 4 4 response 0.74323712 0.04426703 0.64753485 0.8201705 29.0
## 5 5 response 0.92058988 0.03017508 0.83772233 0.9630094 41.0
As the apache
score goes up, the predicted probability of death also goes up. This is what we would expect to see given the positive slope on apache
in the model apache.1
. (Again, if your model had the columns of the response the other way around, you were predicting the probability of survival, and your predictions should then show it going down rather than up to go with the negative slope you would then have had.)
\(\blacksquare\)
- Make a plot of predicted death probability against
apache
score (joined by lines) with, also on the plot, the observed proportion of deaths within eachapache
score, plotted againstapache
score. Does there seem to be good agreement between observation and prediction? Hint: calculate what you need to from the original dataframe first, save it, then make a plot of the predictions, and then to the plot add the appropriate thing from the dataframe you just saved.
Solution
All right, following the hint, let’s start with the original dataframe, called icu
. The plot of the predictions is going to show predicted probabilities of death, so from the data we need to find the observed proportions of death at each apache
score. Our dataframe has one row for each apache
score, so this is not too hard. We have the total number of patients at each score, and the number of deaths out of those patients, so the proportion that died is the second of those divided by the first one:
%>% mutate(proportion = deaths / patients) -> d
icu d
## # A tibble: 38 × 4
## apache patients deaths proportion
## <dbl> <dbl> <dbl> <dbl>
## 1 0 1 0 0
## 2 2 1 0 0
## 3 3 4 1 0.25
## 4 4 11 0 0
## 5 5 9 3 0.333
## 6 6 14 3 0.214
## 7 7 12 4 0.333
## 8 8 22 5 0.227
## 9 9 33 3 0.0909
## 10 10 19 6 0.316
## # … with 28 more rows
That looks all right.
To plot predicted values, your first port of call should be plot_cap
, since the job of this function is to make a lot of predictions and construct a nice plot of them (and so is better than making a plot out of the predictions you made above):
plot_cap(apache.1, condition = "apache")
This again shows that the predicted probability of death goes up sharply with apache
score. Let’s add the data to this graph. The trick is to remember how to add points to a graph you already have, when the points come from a different data set:
plot_cap(apache.1, condition = "apache") +
geom_point(data = d, aes(x = apache, y = proportion), inherit.aes = FALSE)
That does a reasonably good job, but we can do better. We observe now that most of the points are reasonably close to the curve, except for that one in the bottom right, a patient who had the highest apache
score of all, but who happened to survive. The observed proportions that are 1 or 0 over on the right might have been based on only one patient, and the others might have been based on more. So it would be better to encode how many patients each point is based on, for example by the size of the point (you may be able to think of other ways like colour, and you can experiment to see what gets the message across best in your opinion).
A reminder that the inherit.aes
is to tell ggplot
not to take anything for the geom_point
from the plot_cap
graph (which has additional things that it would be a pain to override).
All right, we can make the points different sizes according to the number of patients
they are based on, thus:
plot_cap(apache.1, condition = "apache") +
geom_point(data = d, aes(x = apache, y = proportion, size = patients),
inherit.aes = FALSE)
This shows clearly that the points on the left and right of the graph are based on very few patients, maybe only one each time. But most of the scores between about 10 and 20 were represented by more patients, maybe up to 30. The bigger circles seem to follow the trend of the predictions pretty well, which is what we were hoping for; the proportions based on small numbers of patients might be further away from the predictions, and that’s all right.
Extra: you can see that we want the proportions to be based on reasonable numbers of patients, but the other end of the issue is that we don’t want to combine patients with very different apache
scores, because then you wouldn’t get much of a picture of how well the data follow the trend. This is very much the same idea as choosing the number of bins on a histogram; if you have too many bins, each one will contain very few observations and you won’t see the pattern very clearly, but if you have too few bins, you’ll have lots of observations in each bin but you still won’t see the shape.
With that in mind, perhaps we can try binning the observations in our data here. Let’s see what that does. The starting point is to redefine what I called d
before which had the proportion in it:
<- seq(-1, 45, 4)
break_points break_points
## [1] -1 3 7 11 15 19 23 27 31 35 39 43
%>% mutate(bins = cut(apache, breaks = break_points)) icu
## # A tibble: 38 × 4
## apache patients deaths bins
## <dbl> <dbl> <dbl> <fct>
## 1 0 1 0 (-1,3]
## 2 2 1 0 (-1,3]
## 3 3 4 1 (-1,3]
## 4 4 11 0 (3,7]
## 5 5 9 3 (3,7]
## 6 6 14 3 (3,7]
## 7 7 12 4 (3,7]
## 8 8 22 5 (7,11]
## 9 9 33 3 (7,11]
## 10 10 19 6 (7,11]
## # … with 28 more rows
My bins are 4 apache
points wide. I have no idea yet whether that’s any good, but I want to get the procedure down so that I can come back and change that later if I want to.
What I am using to make the bins is a base R function called cut
. This makes categorical bins out of a quantitative variable (which is normally a bad idea, but we have our reasons here). It takes two inputs: a vector of numbers to bin, and the points that divide one bin from the next, which I defined into a vector called break_points
first. The function cut
defines bins as “half-open”, meaning that the top end is included but the bottom end is excluded. (Note for example which bin an apache
score of 3 goes into.)
Now we want to total up the patients and deaths within each bin, and, for reasons you’ll see later, we want to know where the middle of each bin is (for which I will use the median of the apache
scores within that bin):
%>% mutate(bins = cut(apache, breaks = break_points)) %>%
icu group_by(bins) %>%
summarize(patients = sum(patients), deaths = sum(deaths), apache = median(apache))
## # A tibble: 11 × 4
## bins patients deaths apache
## <fct> <dbl> <dbl> <dbl>
## 1 (-1,3] 6 1 2
## 2 (3,7] 46 10 5.5
## 3 (7,11] 105 19 9.5
## 4 (11,15] 92 32 13.5
## 5 (15,19] 85 36 17.5
## 6 (19,23] 57 34 21.5
## 7 (23,27] 36 23 25.5
## 8 (27,31] 18 12 29.5
## 9 (31,35] 6 6 33.5
## 10 (35,39] 2 2 36.5
## 11 (39,43] 1 0 41
We have redefined the names: patients
and deaths
are now the totals of patients and deaths within each bin, and apache
is now something like the middle apache
score in each bin.
Now we can work out the proportion of patients that died within each bin, and save that:
%>% mutate(bins = cut(apache, breaks = break_points)) %>%
icu group_by(bins) %>%
summarize(patients = sum(patients), deaths = sum(deaths), apache = median(apache)) %>%
mutate(proportion = deaths / patients) -> d
d
## # A tibble: 11 × 5
## bins patients deaths apache proportion
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 (-1,3] 6 1 2 0.167
## 2 (3,7] 46 10 5.5 0.217
## 3 (7,11] 105 19 9.5 0.181
## 4 (11,15] 92 32 13.5 0.348
## 5 (15,19] 85 36 17.5 0.424
## 6 (19,23] 57 34 21.5 0.596
## 7 (23,27] 36 23 25.5 0.639
## 8 (27,31] 18 12 29.5 0.667
## 9 (31,35] 6 6 33.5 1
## 10 (35,39] 2 2 36.5 1
## 11 (39,43] 1 0 41 0
Now we redo our graph, but using the proportions in here as the observed proportions of deaths in the data. You might now realize the reason for those values in the apache
column: on the graph, the proportion
values will be our \(y\) coordinates, but we needed something to be the \(x\) coordinates.
Actually redrawing our graph is in fact exactly the same code as before; the thing that has changed is our definition of d
in it:
plot_cap(apache.1, condition = "apache") +
geom_point(data = d, aes(x = apache, y = proportion, size = patients),
inherit.aes = FALSE)
And now we see, with more patients per bin, that the agreement between data and prediction is very good apart from the very small bins on the right.
One modification that you might like to pursue is to have variable-width bins: wider bins at the extremes and narrower bins in the middle, so that each bin has about the same number of patients. You could do this by modifying what I called break_points
so that the numbers in it are more spread out at the extremes and closer together in the middle.
Extra extra: if you’re a football fan, you could imagine doing a very similar analysis by modelling the probability of successfully kicking a field goal as it depends on the distance it is attempted from. This would probably work pretty well for the NFL and (no doubt) for the CFL as well. In both these codes of football, field goals are always attempted from between two side-to-side marks on the field called “hash marks”. In the NFL, the hash marks are close together, so field goals are attempted from more or less directly in front of the goalposts. In the CFL, the hash marks are further apart, so a field goal might be attempted from off to the side of the goalposts, and then the probability of success might also depend on how far off to the side you are.
If you happen also to be a rugby fan, you’ll know that kicks at goal might need to be attempted from anywhere on the field, even out near the sidelines, and in that case, the chance of kicking a goal definitely depends on where on the field you are kicking it from as well as how far out it is.
Extra extra extra: football (NFL and CFL) have an effective formality called the “point after”: after a team scores a touchdown, the kicker kicks at goal from right in front of the posts, which is almost never missed. In rugby, the equivalent to a touchdown is called a “try”, and to score a try the player must literally touch the ball down. (Thus, to a rugby fan, the NFL and CFL “touchdowns” are absurd misnomers.) After a try is scored, there is a kick at goal (called a “conversion”) which is not taken from right in front of the posts like in the NFL, but from the side-to-side position on the field where the try was scored. Thus, if a try is scored near the sidelines, as it often is, the conversion has to be kicked from near the sideline as well. The kicker can go as far back down the field as they like, but as you might imagine, a sideline conversion is still very difficult to kick.
Extra-to-the-fourth: the reason for those names in rugby. In the old days, if a team scored a try and kicked the conversion as well, it was called a “goal”, and as a rugby player, you wanted your team to score a goal. Scoring a try allowed a team to “try” to “convert” it into a “goal” by kicking the conversion, hence the names. Because the NFL and CFL have their roots in rugby, they have the “point after” a touchdown as well (with the option now of allowing teams to try for two points after a touchdown as well).
\(\blacksquare\)
26.20 Go away and don’t come back!
When a person has a heart attack and survives it, the major concern of health professionals is to prevent the person having a second heart attack. Two factors that are believed to be important are anger and anxiety; if a heart attack survivor tends to be angry or anxious, they are believed to put themselves at increased risk of a second heart attack.
Twenty heart attack survivors took part in a study. Each one was given
a test designed to assess their anxiety (a higher score on the test
indicates higher anxiety), and some of the survivors took an anger
management course.
The data are in
link; y
and
n
denote “yes” and “no” respectively. The columns denote
(i) whether or not the person had a second heart attack, (ii) whether
or not the person took the anger management class, (iii) the anxiety
score.
- Read in and display the data.
Solution
<- "http://ritsokiguess.site/datafiles/ha2.txt"
my_url <- read_delim(my_url, " ") ha
## Rows: 20 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: " "
## chr (2): second, anger
## dbl (1): anxiety
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ha
## # A tibble: 20 × 3
## second anger anxiety
## <chr> <chr> <dbl>
## 1 y y 70
## 2 y y 80
## 3 y y 50
## 4 y n 60
## 5 y n 40
## 6 y n 65
## 7 y n 75
## 8 y n 80
## 9 y n 70
## 10 y n 60
## 11 n y 65
## 12 n y 50
## 13 n y 45
## 14 n y 35
## 15 n y 40
## 16 n y 50
## 17 n n 55
## 18 n n 45
## 19 n n 50
## 20 n n 60
The anxiety scores are numbers; the other two variables are “yes” and “no”, which makes perfect sense.
\(\blacksquare\)
- * Fit a logistic regression predicting whether or not a heart attack survivor has a second heart attack, as it depends on anxiety score and whether or not the person took the anger management class. Display the results.
Solution
This:
.1 <- glm(factor(second) ~ anger + anxiety, family = "binomial", data = ha)
hasummary(ha.1)
##
## Call:
## glm(formula = factor(second) ~ anger + anxiety, family = "binomial",
## data = ha)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.52106 -0.68746 0.00424 0.70625 1.88960
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.36347 3.21362 -1.980 0.0477 *
## angery -1.02411 1.17101 -0.875 0.3818
## anxiety 0.11904 0.05497 2.165 0.0304 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.726 on 19 degrees of freedom
## Residual deviance: 18.820 on 17 degrees of freedom
## AIC: 24.82
##
## Number of Fisher Scoring iterations: 4
There’s no need to worry about two-column responses here, because each
row of the data frame refers to only one person, and the response
variable second
is already a categorical variable with two
categories (that needs to be turned into a factor
for
glm
).
\(\blacksquare\)
- In the previous part, how can you tell that you were predicting the probability of having a second heart attack (as opposed to the probability of not having one)?
Solution
The levels of a factor are taken in alphabetical order, with
n
as the baseline, so we are predicting the probability
of the second one y
.
\(\blacksquare\)
- *
For the two possible values
y
andn
ofanger
and the anxiety scores 40, 50 and 60, make a data frame containing all six combinations, and use it to obtain predicted probabilities of a second heart attack. Display your predicted probabilities side by side with what they are predictions for.
Solution
This time, I give you the values I want the predictions for (as opposed to calculating something like quartiles from the data), so we might as well just type them in.
First step, datagrid
:
<- datagrid(model = ha.1, anger = c("y", "n"), anxiety = c(40, 50, 60))
new new
## anger anxiety
## 1: y 40
## 2: y 50
## 3: y 60
## 4: n 40
## 5: n 50
## 6: n 60
Then, the predictions:
predictions(ha.1, newdata = new)
## rowid type predicted std.error conf.low conf.high anger anxiety
## 1 1 response 0.06749763 0.08284878 0.005455627 0.4885217 y 40
## 2 2 response 0.19226955 0.15447608 0.032771639 0.6257931 y 50
## 3 3 response 0.43908386 0.22715044 0.113796084 0.8267515 y 60
## 4 4 response 0.16774789 0.16551787 0.019352747 0.6730543 n 40
## 5 5 response 0.39861864 0.19858379 0.115592812 0.7707225 n 50
## 6 6 response 0.68551307 0.16294573 0.331329889 0.9055632 n 60
Extra: In the help file for predictions
, you will see it done in one step, with code like
predictions(ha.1, newdata = datagrid(anger = c("y", "n"), ...))
This saves constructing new
first, but the code for datagrid
is to my mind usually rather long to fit on one line inside the predictions(...)
.
\(\blacksquare\)
- Use your predictions from the previous part to describe the
effect of changes in
anxiety
andanger
on the probability of a second heart attack.
Solution
The idea is to change one variable while leaving the other fixed. (This is the common refrain of “all else equal”.)
Pick a level of anger
, say n
(it doesn’t matter
which) and look at the effect of anxiety
. The probability
of a second heart attack increases sharply from 0.17 to 0.40 to
0.69. So an increased anxiety score is associated with an
increased probability of second heart attack (all else equal).
To assess the effect of taking the anger management course, pick
an anxiety
value, say 40, and compare the probabilities
for anger
n
and y
. For someone who has
not taken the anger management course, the probability is 0.17,
but for someone who has, it drops all the way to 0.07. (The
pattern is the same at the other anxiety scores.)
Extra: the reason it doesn’t matter what value of the other
variable you look at (as long as you keep it fixed) is that the
model is “additive” with no interaction, so that the effect of
one variable does not depend on the value of the other one. If we
wanted to see whether the effect of anxiety was different
according to whether or not the person had taken the anger
management course, we would add an interaction between
anxiety
and anger
. But we won’t be doing this
kind of thing until we get to analysis of variance, so you don’t
need to worry about it yet.
\(\blacksquare\)
- Are the effects you described in the previous part
consistent with the
summary
output fromglm
that you obtained in (here)? Explain briefly how they are, or are not. (You need an explanation for each ofanxiety
andanger
, and you will probably get confused if you look at the P-values, so don’t.)
Solution
In the previous part, we found that increased anxiety went with an
increased probability of second heart attack. Back in
(here), we found a positive slope of 0.11904 for anxiety,
which also means that a higher anxiety score goes with a higher
probability of second heart attack.
That was not too hard. The other one is a little more work.
anger
is categorical with two categories n
and
y
. The first one, n
, is the baseline, so
angery
shows how y
compares to n
. The
slope \(-1.04211\) is negative, which means that someone who has
taken the anger management course has a lower probability
of a second heart attack than someone who hasn’t (all else
equal). This is the same story that we got from the predictions.
That’s it for that, but I suppose I should talk about those
P-values.
anxiety
is significant, so it definitely has an effect on
the probability of a second heart attack. The pattern is clear
enough even with this small data set. Here’s a visual:
ggplot(ha, aes(x = second, y = anxiety)) + geom_boxplot()
This is the wrong way around in terms of cause and effect, but it shows pretty clearly that people who suffer a second heart attack have a higher level of anxiety than those who don’t.
That’s clear enough, but what about anger
? That is not
significant, but there seems to be a visible effect of anger
on the predicted probabilities of (here): there we saw
that if you had done the anger management course, your probability of
a second heart attack was lower. But that’s only the predicted
probability, and there is also uncertainty about that, probably quite
a lot because we don’t have much data. So if we were to think about
confidence intervals for the predicted probabilities, they would be
wide, and for the two levels of anger
at a fixed
anxiety
they would almost certainly overlap.
Another way of seeing that is a visual, which would be a side-by-side bar chart:
ggplot(ha, aes(x = anger, fill = second)) + geom_bar(position = "dodge")
A small majority of people who took the anger management did not have a second heart attack, while a small minority of those who did not, did.50 But with these small numbers, the difference, even though it points the way we would have guessed, is not large enough to be significant:
with(ha, table(anger, second))
## second
## anger n y
## n 4 7
## y 6 3
This is not nearly far enough from an equal split to be significant. (It’s the same kind of idea as for Mood’s median test in C32/C33.)
\(\blacksquare\)
I get confused about the difference between morals and ethics. [This is a very short description of that difference:] (http://smallbusiness.chron.com/differences-between-ethical-issues-moral-issues-business-48134.html). The basic idea is that morals are part of who you are, derived from religion, philosophy etc. Ethics are how you act in a particular situation: that is, your morals, what you believe, inform your ethics, what you do. That’s why the students had to play the role of an ethics committee, rather than a morals committee; presumably the researcher had good morals, but an ethics committee had to evaluate what he was planning to do, rather than his character as a person.↩
As with many of these acronyms, you get the idea that the acronym came first and they devised some words to fit it.↩
My daughter learned the word pre-empt because we like to play a bridge app on my phone; in the game of bridge, you make a pre-emptive bid when you have no great strength but a lot of cards of one suit, say seven, and it won’t be too bad if that suit is trumps, no matter what your partner has. If you have a weakish hand with a lot of cards in one suit, your opponents are probably going to be able to bid and make something, so you pre-emptively bid first to try and make it difficult for them.↩
This kind of thing is sometimes called an inverse prediction.↩
The usual application of this is to combine a number with a vector. If you try to subtract a length-2 vector from a length-6 vector, R will repeat the shorter one three times and do the subtraction, but if you try to subtract a length-2 vector from a length-7 vector, where you’d have to repeat the shorter one a fractional number of times, R will do it, but you’ll get a warning, because this is probably not what you wanted. Try it and see.↩
This is done in the
geom_point
on the website, which is where I learned about it.↩The test is this way because it’s a generalized linear model rather than a regular regression.↩
See, I told you it was a little confusing.↩
When you have one observation per line, the predictions are of the second of the two levels of the response variable. When you make that two-column response, the predictions are of the probability of being in the first column. That’s what it is. As the young people say, don’t @ me.↩
When you have higher-order terms, you have to keep the lower-order ones as well: higher powers, or interactions (as we see in ANOVA later).↩
I get confused about the difference between morals and ethics. [This is a very short description of that difference:] (http://smallbusiness.chron.com/differences-between-ethical-issues-moral-issues-business-48134.html). The basic idea is that morals are part of who you are, derived from religion, philosophy etc. Ethics are how you act in a particular situation: that is, your morals, what you believe, inform your ethics, what you do. That’s why the students had to play the role of an ethics committee, rather than a morals committee; presumably the researcher had good morals, but an ethics committee had to evaluate what he was planning to do, rather than his character as a person.↩
If you took STAB23, you’ll have used PSPP, which is a free version of SPSS.↩
If it is for you, go right ahead, but for me it wasn’t.↩
As with many of these acronyms, you get the idea that the acronym came first and they devised some words to fit it.↩
Read that carefully.↩


The Sweet Hereafter
(Canada, 1997)
Seen on 24 December 1998 with Linda for $8.75 at the Lincoln Plaza Cinema
Perhaps I am turning into a stone, but I was unmoved by *The Sweet Hereafter*, although it is a beautifully filmed movie. Other reviewers have described it as elegaic. While it does remind you of an elegy, how many elegies run feature-length?
*The Sweet Hereafter* examines the lives of parents just after a fatal accident in which a bus plunges through a guard rail, careens down an embankment, and falls through the frozen surface of a lake in British Columbia.
The story that drives the movie is the arrival of an attorney, Mitchell Stephens (Ian Holm) who wants to push some of the grieving parents into pursuing compensation, that perhaps somehow the senseless tragedy can be blamed on someone or something; a careless school bus driver, or a faulty bolt, or something. As if monetary compensation will make up for the loss of fourteen children. This sort of lawsuit will gain him one-third of the netted compensation. Also, it means he has to find the "right" sort of pa rents to testify.
There are also two parallels to the main storyline. First, there is the irony of Stephens appealing to parents when his own daughter, a drug addict, is out of control and he can no longer do anything for her. Or, as he put it, his love for her has become and anger something like "steaming piss." Also, there is the parallel of the children's story of the Pied Piper of Hamelin. This is one of the strongest points of the movie. A lot of the parents in *The Sweet Hereafter* are hardly worthy of their children, and there is a strong sense of not just sadness but punishment as well.
Director Atom Egoyan does a good job of showing people after tragedy. There is a strong sense of the Time Before and the Time After. Unfortunately, from the beginning of the movie, it is difficult to follow the timing of the story and the many flashbacks, as there is a lot of jumping around. Holm is shown having a chance meeting on a plane--a friend of his daughter--to whom he talks of his daughter's troubles. Where that fits into the timing of the rest of the story is very difficult.
I think that because the story is told subtly, and because the characters seem so low-key and almost stoic, it is hard to really empathize with them. There is also the taboo of incest, a topic that should have evoked more emotion, but Egoyan does a lot of filmus interruptus, and the constant jumping around of flashbacks hurts the full development of investing in the story and the characters. You are on the edge of tears, and then suddenly, you are pulled out of the moment. Even when the characters are fighting, as in the scene where Billy Ansell (Bruce Greenwood) and Sam Burnell (Tom McCamus) disagree over having a group lawsuit, there is no intensity.
A town losing 14 children in a tragedy undoubtedly changes that town forever, and Egoyan shows the emptiness and sterility very well. It's just a shame that it was not as moving as it could have been. *The Sweet Hereafter* left me with the feeling of looking at the carcass of a wrecked bus, but not knowing anything about the people who perished within or the people who loved them. Feeling sorry for them is as difficult as feeling sorry for the people of Hamelin. There is more after than before to really build the sort of sympathy a movie like this surely wants to achieve.
Written by Atom Egoyan, based on a novel by Russell Banks
More movie reviews by Seth Bookey, with graphics, can be found at http://www.geocities.com/Athens/2679/kino.html
The review above was posted to the
rec.arts.movies.reviews newsgroup (de.rec.film.kritiken for German reviews).
The Internet Movie Database accepts no responsibility for the contents of the
review and has no editorial control. Unless stated otherwise, the copyright
belongs to the author.
Please direct comments/criticisms of the review to relevant newsgroups.
Broken URLs inthe reviews are the responsibility of the author.
The formatting of the review is likely to differ from the original due
to ASCII to HTML conversion.
Related links: index of all rec.arts.movies.reviews reviews


W
idth
L
ength 12" or 14
-
1/2"
31" $96
36" $98
39" $105
43" $107
48" $114
52" $118
55" $123
60" $129
64" $145
67" $152
72" $158
75" $165
80" $167
raised pan
e
l shutters
INDEX
Custom Sizes Height must be special quoted


My name is Wyatt Salmons. I'm 28 years old, and from Albuquerque, New Mexico. I have 2 dogs and 3 cats (but as for the cats I wish I didn't lol). I have lived in Belen, NM for about 2 years with my fiance, Jessica. We love it! Its quiet, cheap, and slow.


Generates a calendar along the lines of this one (with month and day names in the configured locale):
< January 2003 >
Mo Tu We Th Fr Sa Su
1 2 3 4 5
6 7 8 9 10 11 12
13 14 15 16 17 18 19
20 21 22 23 24 25 26
27 28 29 30 31
It walks through all your entries and marks the dates that have entries so you can click on the date and see entries for that date.
This plugin comes with Pyblosxom. To install, do the following:
You can set the start of the week using the calendar_firstweekday configuration setting, for example:
py['calendar_firstweekday'] = 0
will make the week start on Monday (day ‘0’), instead of Sunday (day ‘6’).
Pycalendar is locale-aware. If you set the locale config property, then month and day names will be displayed according to your locale.
It uses the following CSS classes:
Plugin is distributed under license: Public domain


Hi! I'm Paras Vora.
I'm an engineering student who, for some reason, wants to be a physician.
I care about making the world a better place. To do that, I develop and apply technological solutions to pre-existing problems.
Check out my Résumé View My Portfolio Get in TouchWhat I'm about
I'm a 20-year old student from Kentucky, and spend my time learning, researching potential cancer treatments, helping fundraise for the American Cancer Society, playing tennis, working on small web projects, and amateur photography.
Combined Passions
As a biomedical engineer, computer science major, and pre-medical student, I have a range of experiences to draw upon, and want to use the skills I have developed in order to help people in a meaningful way.
Easy to work with
I'm capable of working with, and leading, groups of all sizes. I make sure that my documentation is detailed and specific enough to help anyone get started. I'm open to feedback and criticism
I love innovation and admire creativity. As a future physician, I want to take my background and apply it to a field that is in constant need of innovation, with a little bit of creativity. My passion for problem-solving drives me to finding lasting solutions to current problems.
Portfolio
A team project, using Java and SQL, and the Eye Tribe eye-tracker to evaluate reading performance. Competed in a poster competition.
As part of my rapid prototyping course's creative programming assignment, and my involvement with Relay For Life, I revamped our check-in process.
Assisted a research project investigating the effects of novel anti-cancer drugs on bone homeostasis. Conducted bone histomorphometric analysis.
Interested in working with me?
Check out my Résumé Get in TouchI hope my skillset can help you achieve your goals in making a difference.


BLM Intel: Real time protest information and resources.
Home
More Info
Want to help out or learn more? Check out the links below for more information!
Send an email to demand justice for George Floyd.
Lots of information/resources
George Floyd Petition
List of peitions you can sign
Another huge list of resources
Find protests in your area
Another great protest tracker


Use choose sig file to pick a sig file
Use choose data file to pick a test data file
Use Test Signature to test the sig against your test data


This is my first webpage!
Isn't it nice?
here are some things I'm going to learn about coding:
Check out Epicodus so you can make your own!
By the way, here are some of my favorite things.
I also have some favorite places if you've got some extra time to browse around.


| BII2-T10-R001 |
BII2 |
|
cbc:CustomizationID |
| SI-V11-PO-R200 |
SI |
This XML instance is NOT tagged as an SI-UBL order; please check the CustomizationID value |
contains(normalize-space(.), 'urn:www.cenbii.eu:transaction:biitrns001:ver2.0:extended:urn:www.peppol.eu:bis:peppol3a:ver2.0:extended:urn:www.simplerinvoicing.org:si:si-ubl:ver1.2') |
| SI-PO-R000 |
SI |
This is NOT a UBL 2.x Order, validation cannot continue |
ubl:Order |
| SI-V12-PO-R201 |
SI |
This XML instance is NOT tagged as an SI-UBL order |
contains(., 'urn:www.simplerinvoicing.org:si:si-ubl:ver1.2') |
| SI-V12-PO-R001 |
SI |
Each order MUST be specified as a UBL 2.1 document |
cbc:UBLVersionID='2.1' |
| SI-V12-CODE-R003 |
SI |
Country codes in an order MUST be coded using ISO code list 3166-1 |
Long XPath used, check the schematron source |
| BII2-T01-R007 |
BII2 |
Allowances and charges MUST have a reason |
(cbc:AllowanceChargeReason) |
| BII2-T01-R021 |
BII2 |
An order MUST have the buyer party name or a buyer party identifier |
(cac:Party/cac:PartyName/cbc:Name) or (cac:Party/cac:PartyIdentification/cbc:ID) |
| BII2-T01-R012 |
BII2 |
Expected total amount for payment MUST NOT be negative, if expected total amount for payment is provided |
(cbc:PayableAmount) >= 0 |
| BII2-T01-R013 |
BII2 |
Expected total sum of line amounts MUST NOT be negative, if expected total sum of line amounts is provided |
(cbc:LineExtensionAmount) >= 0 |
| BII2-T01-R014 |
BII2 |
Expected total sum of line amounts MUST equal the sum of the order line amounts at order line level, if expected total sum
of line amounts is provided
|
(xs:decimal(cbc:LineExtensionAmount)) = (round(xs:decimal(sum(/ubl:Order/cac:OrderLine/cac:LineItem/cbc:LineExtensionAmount))
* 10 * 10) div 100)
|
| BII2-T01-R015 |
BII2 |
Expected total sum of allowance at document level MUST be equal to the sum of allowance amounts at document level, if expected
total sum of allowance at document level is provided
|
(cbc:AllowanceTotalAmount and (xs:decimal(cbc:AllowanceTotalAmount)) = (round(xs:decimal(sum(/ubl:Order/cac:AllowanceCharge[cbc:ChargeIndicator="false"]/cbc:Amount))
* 10 * 10) div 100)) or not(cbc:AllowanceTotalAmount)
|
| BII2-T01-R016 |
BII2 |
Expected total sum of charges at document level MUST be equal to the sum of charges at document level, if expected total sum
of charges at document level is provided
|
(cbc:ChargeTotalAmount and (xs:decimal(cbc:ChargeTotalAmount)) = (round(xs:decimal(sum(/ubl:Order/cac:AllowanceCharge[cbc:ChargeIndicator="true"]/cbc:Amount))
* 10 * 10) div 100)) or not(cbc:ChargeTotalAmount)
|
| BII2-T01-R017 |
BII2 |
Expected total amount for payment MUST be equal to the sum of line amounts minus sum of allowances at document level plus
sum of charges at document level and VAT total amount, if expected total amount for payment is provided
|
Long XPath used, check the schematron source |
| BII2-T01-R019 |
BII2 |
Each item property MUST have a data name, if item property is provided |
(cbc:Name) |
| BII2-T01-R020 |
BII2 |
Each item property MUST have a data value, if item property is provided |
(cbc:Value) |
| BII2-T01-R001 |
BII2 |
An order MUST have a customization identifier |
(cbc:CustomizationID) |
| BII2-T01-R002 |
BII2 |
An order MUST have a profile identifier |
(cbc:ProfileID) |
| BII2-T01-R004 |
BII2 |
An order MUST have a document issue date |
(cbc:IssueDate) |
| BII2-T01-R006 |
BII2 |
An order MUST have a document identifier |
(cbc:ID) |
| BII2-T01-R009 |
BII2 |
An order MUST be stated in a single currency |
(cbc:DocumentCurrencyCode) |
| BII2-T01-R003 |
BII2 |
Each order line MUST have a document line identifier that is unique within the order |
(cac:LineItem/cbc:ID) |
| BII2-T01-R010 |
BII2 |
Each order line ordered quantity MUST not be negative |
(cac:LineItem/cbc:Quantity) >= 0 |
| BII2-T01-R011 |
BII2 |
Each order line item net price MUST not be negative |
((cac:LineItem/cac:Price/cbc:PriceAmount) >= 0) or (not(cac:LineItem/cac:Price)) |
| BII2-T01-R030 |
BII2 |
Each order line ordered quantity MUST have an associated unit of measure |
(cac:LineItem/cbc:Quantity/@unitCode) |
| BII2-T01-R031 |
BII2 |
Each order line MUST have an item identifier and/or an item name |
(cac:LineItem/cac:Item/cbc:Name) or (cac:LineItem/cac:Item/cac:StandardItemIdentification/cbc:ID) or (cac:LineItem/cac:Item/cac:SellersItemIdentification/cbc:ID) |
| BII2-T01-R022 |
BII2 |
An order MUST have the seller party name or a seller party identifier |
(cac:Party/cac:PartyName/cbc:Name) or (cac:Party/cac:PartyIdentification/cbc:ID) |
| EUGEN-T01-R015 |
EUGEN |
A country
identification code MUST have a list identifier attribute
“ISO3166-1:Alpha2”
|
@listID = 'ISO3166-1:Alpha2' |
| EUGEN-T01-R014 |
EUGEN |
A document
currency code MUST have a list identifier attribute “ISO4217”
|
@listID='ISO4217' |
| EUGEN-T01-R011 |
EUGEN |
An endpoint
identifier MUST have a scheme identifier attribute
|
@schemeID |
| EUGEN-T01-R013 |
EUGEN |
An order type
code MUST have a list identifier attribute “UNCL1001”
|
@listID='UNCL1001' or @listID='UNCLD1001' |
| EUGEN-T01-R012 |
EUGEN |
A party
identifier MUST have a scheme identifier attribute
|
@schemeID |
| EUGEN-T01-R017 |
EUGEN |
A tax category
identifier MUST have a scheme identifier attribute “UNCL5305”
|
@schemeID = 'UNCL5305' |
| EUGEN-T01-R016 |
EUGEN |
A unit code
attribute MUST have a unit code list identifier attribute “UNECERec20”
|
not(attribute::unitCode) or (attribute::unitCode and attribute::unitCodeListID = 'UNECERec20') |
| OP-T01-R001 |
OP |
OrderTypeCode in an order MUST be from UN/ECE 1001 code list PEPPOL Subset |
Long XPath used, check the schematron source |
| OP-T01-R002 |
OP |
DocumentCurrencyCode MUST be coded using ISO code list 4217 |
Long XPath used, check the schematron source |
| OP-T01-R003 |
OP |
A Binary Object MIME code attribute MUST be coded using MIME Code Type version 2008 code list |
Long XPath used, check the schematron source |
| OP-T01-R004 |
OP |
An Enpoint Identifier Scheme identifier MUST be from the list of PEPPOL Party Identifiers described in the "PEPPOL Policy
for using Identifiers".
|
Long XPath used, check the schematron source |
| OP-T01-R005 |
OP |
An Party Identifier Scheme identifier MUST be from the list of PEPPOL Party Identifiers described in the "PEPPOL Policy for
using Identifiers".
|
Long XPath used, check the schematron source |
| OP-T01-R006 |
OP |
Country codes in an order MUST be coded using ISO code list 3166-1 |
Long XPath used, check the schematron source |
| OP-T01-R007 |
OP |
currencyID MUST be coded using ISO code list 4217 |
Long XPath used, check the schematron source |
| OP-T01-R008 |
OP |
A tax category identifier MUST be coded using UN/ECE 5305 BII2 Subset |
( ( not(contains(normalize-space(.),' ')) and contains( ' AE E S Z AA H ',concat(' ',normalize-space(.),' ') ) ) ) |
| OP-T01-R009 |
OP |
Unit code MUST be coded according to the UN/ECE Recommendation 20 (2009) |
Long XPath used, check the schematron source |


Add expenses
Add income
Description:
Cost
Category
School costs
Entertainment
Food
Housing
Transportation
Loans/Payments
Other
Add
Description:
Amount
Category
Wages
Financial aid
From savings
Assistance from parents
Other
Add


Alerts are added and configured in the Alert Tab of any dashboard graph panel, letting you build and visualize an alert using existing queries. To persist your alert rule changes remember to save the dashboard.


Get Started
We provide open source APIs to get hospital bed status in Banglore region. We will extend our APIs to different cities and add further in the coming weeks to provide useful data to the developer community.
All requests should be directed to sevaprod.sagisu.com (Auth tokens not required)
{{url}} = https://sevaprod.sagisu.com/v1
/bedStatus
GET : {{url}}/bedStatus
Provides bed count for hospitals listed on the website. API sorts in descending order of total beds available in each hospital.
Note : Query param can be used to sort on ICU count or ICU ventilator count
Query Params
sortKey : icu or icuVentilator
sortOrder : 1(ascending) or -1(descending)
Sample request : https://sevaprod.sagisu.com/v1/bedStatus
Sample response :
{
"data": [
{
"_id": "6092345687e565515850cb7d",
"name": "ESIMC PGIMSR MODEL HOSPITAL",
"address": " 620/A/1, 12th Main Rd, 2nd Block, Rajajinagar, Bengaluru, Karnataka 560010, India",
"mobile": "<PHONE_NUMBER>30",
"general": 11,
"hdu": 0,
"icu": 0,
"icuVentilator": 0,
"total": 11
}
]
}
Response fields
| Key | Description |
|---|---|
| _id | Unique ID of the hospital |
| name | Name of the hospital |
| address | Address of the hospital |
| mobile | Mobile number of the hospital |
| general | Number of beds available in general ward |
| hdu | Number of beds available in HDU ward |
| icu | Number of beds available in ICU |
| icuVentilator | Number of beds available in ICU ventilator |
| total | Total number of beds available in the hospital across all the wards |
POST : {{url}}/bedStatus
API to add new hospital details to the system or to update existing hospital details.Accepts array of hospital details to be added or corrected in the system.
Sample request : https://sevaprod.sagisu.com/v1/bedStatus
Body params : array of json objects
[
{
"name": "General Hospital Jayanagar",
"address": "32nd Cross, 4th B Block, beside Rajiv Gandhi Health Institute, Tilak Nagar, Jayanagar, Bengaluru, Karnataka 560041, India",
"mobile": "918022445711",
"general": 6,
"hdu": 0,
"icu": 0,
"icuVentilator": 0,
"total": 6
},
{
"name": "General Hospital Anekal",
"address": "Weaver Colony, Anekal, Karnataka 562106, India",
"general": 18,
"hdu": 0,
"icu": 0,
"icuVentilator": 0,
"total": 18
}
]
Request fields
| Key | Validation | Description |
|---|---|---|
| name | required,string | Name of the hospital |
| address | string | Address of the hospital |
| mobile | string | Mobile number of the hospital |
| general | Number | Number of beds available in general ward |
| hdu | Number | Number of beds available in HDU ward |
| icu | Number | Number of beds available in ICU |
| icuVentilator | Number | Number of beds available in ICU ventilator |
| total | Number | Total number of beds available in the hospital across all the wards |
/reserve
GET : {{url}}/reserve
Gives the list of bookings made by the user
Query Params
userId : id of the user recieved during login
Sample request : https://sevaprod.sagisu.com/v1/reserve?userId=60b90516499210371968513f
Sample response :
{
"data":[
{
"id":"60f8fdabb2c8d909bca63dae",
"reserve":{
"id":"60f8fdabb2c8d909bca63db0",
"type":"ICU"
},
"status":"PENDING",
"patientName":"Test User",
"age":"20",
"location":"User location",
"onOxygenCylinder":false,
"hospitalsVisited":[A,B,C],
"testedForCovid":true,
"covidResult":"Negative",
"hospitalPreference":"Private",
"attenderMobile":"XXXXXXXXXX",
"relationship":"wife",
"hospital":"60ee74d2973f0b38081b72f7",
"user":"60b90516499210371968513f"
}
]
}
Response fields
| Key | Description |
|---|---|
| id | Unique ID of the reservation request |
| reserve.type | Type of the hospital bed booked |
| status | Status of reservation request(PENDING,ACCEPTED,REJECTED). |
| hospital | ID of the hospital booked |
| user | ID of the user |
| Other patient related informations |
POST : {{url}}/reserve
API to make hospital bed bookings.
Sample request : https://sevaprod.sagisu.com/v1/reserve
Body params : json of booking information
{
"reserve":{
"type":"ICU"
},
"patientName":"Test User",
"age":"20",
"location":"User location",
"onOxygenCylinder":false,
"hospitalsVisited":[A,B,C],
"testedForCovid":true,
"covidResult":"Negative",
"hospitalPreference":"Private",
"attenderMobile":"XXXXXXXXXX",
"relationship":"wife",
"srfid":"900",
"buNumber" : "900",
"user":"60b90516499210371968513f"
}
Request fields
| Key | Validation | Description |
|---|---|---|
| reserve.type | required,string | Type of the bed requested |
| patientName | required,string | Name of the patient |
| age | required,string | Age of the patient |
| location | required,string | Address of the patient |
| onOxygenCylinder | required,boolean | Is pateint on oxygen cylinder or not |
| hospitalsVisited | string array | List of all the hospitals visited before booking here |
| testedForCovid | required,boolean | Is covid test is done or not? |
| covidResult | required,string | Covid report details. POSITIVE,NEGATIVE,AWAITING are the values to this field |
| hospitalPreference | required,string | Preferences for the hospital. GOVERNMENT,PRIVATE or ANY should be the values to this field |
| attenderMobile | required,string | Mobile number of the attender to which the booking details will be sent. |
| relationship | required,string | Relaionship of the attender with the patient |
| srfid | string | SRFID given during covid test |
| buNumber | string | BU Number given during covid test |
| user | string | ID of the user, you will get it during login |


Datasets used in "Statistical Methods for the Social Sciences" (SMSS) by Alan Agresti and Barbara Finlay.
| Version: | 1.0-2 |
| Suggests: | devtools |
| Published: | 2015-10-31 |
| Author: | Jeffrey B. Arnold [aut, cre], Alan Agresti [cph], Barbara Finlay [cph] |
| Maintainer: | Jeffrey B. Arnold <jeffrey.arnold at gmail.com> |
| License: | GPL-3 |
| NeedsCompilation: | no |
| Materials: | README NEWS |
| CRAN checks: | smss results |
| Reference manual: | smss.pdf |
| Package source: | smss_1.0-2.tar.gz |
| Windows binaries: | r-devel: smss_1.0-2.zip, r-release: smss_1.0-2.zip, r-oldrel: smss_1.0-2.zip |
| OS X binaries: | r-release: smss_1.0-2.tgz, r-oldrel: smss_1.0-2.tgz |
| Old sources: | smss archive |
Please use the canonical form https://CRAN.R-project.org/package=smss to link to this page.


"That awful little man on your shoulder"
June 5, 2010
One of my very favourite Brooklyn-based producers, Matt Shadetek, wrote this post about creativity and the fears our society, starting by ourselves, imposes to be truly creative. I recommend a full read of it because it's spotted with very compelling thoughts and honest acknowledgements about his experience in releasing personal music. This one, particularly, made me feel identified:
(...) The problem for me is that fear, especially in releasing an instrumental which I produced alone, with no collaborating producers or vocalists, that in doing so there is nothing to hide behind and people will not just judge my music but they will judge me. (...)
I'm nothing like a musician or an artist, but I've always liked the idea of being a PhD student as something not that far: art or not, research is a creative activity and, whether or not the results are as enoyable or useful as those from Shadetek, for instance, there are several similarities in the processes of creating songs and academic papers.
I am about to hand in the first final draft of my dissertation, and these days the following question starts popping up in conversations more often: "Oh, that's cool you are finishing your dissertation, what is it about?" I've always had a very hard time answering it. And it's not because the dissertation topic (cities and their spatial structure) is hard to explain to a non-expert on the field; although there might be some of it, I think most of the trouble relates to Shadetek's point: it's the fear of people thinking what I do is not relevant, useful or creative enough, of people going beyond and not judging my work but judging me that usually takes me to the easy way of "I don't wanna bore you now, I'll tell you other day if you're interested".   They’re all going to laugh at you'”.
For those not that interested so much in my PhD's suffering, I've also got a paragraph for you. I think in the end, whether it's a song, a dissertation or your love for 1950's icelandic folk music, what makes them scaring is that they are all things that relate to our very inner self and, by exposing them to others we feel somehow naked, afraid of being neglected because we are in a way different or particular, and it's in our human condition to seek for belonging and social acceptance. Once I heard that "shyness is the desire to please and the fear not to success on it". At least for myself, that's the little man on my shoulder.


Scott Doyle /<EMAIL_ADDRESS>Good data (descriptive features, low noise, high $N$) means that any classifier will do a decent job.
However, some are more suited to different situations than others.
Posterior probability:
$ P(\omega_{j}|\mathbf{x}) = \frac{p(\mathbf{x}|\omega_{j})P(\omega_{j})}{p(\mathbf{x})}$
$ p(\mathbf{x}) = \sum_{j=1}^{c}p(\mathbf{x}|\omega_{j})P(\omega_{j}) $
In bagging, training subsets of $\mathcal{D}$ are created by drawing $n^{\prime}<n$ samples with replacement.
Each component classifier casts a “vote” for the classification of a sample point, and the final classification is the result of this vote.
Each component classifier can be the same type, although the learned parameter values may vary (since they are created using different training sets).
This increases the stability of the final classifier by averaging over the differences incurred by using different training sets.
In boosting, the goal is to increase the accuracy of the final classifier.
In this case, we use component classifiers, each of which is the “most informative” out of all possible component classifiers.
So we are less concerned about reducing “noise” or increasing stability; here, we want to combine different aspects of the dataset (represented by classifiers) to reach the right decision.
Two-category problem with three component classifiers: $C_{1}, C_{2}, C_{3}$
Each $C_{i}$ is trained on training sets $\mathcal{D}_{1}, \mathcal{D}_{2}, \mathcal{D}_{3}$.
We select $\mathcal{D}_{2}$ such that half of its samples are misclassified by $C_{1}$.
$\mathcal{D}_{3}$ contains samples where $C_{1}$ and $C_{2}$ disagree.
How should we select samples for $\mathcal{D}_{1}$?
Ideally, $n_{1}\simeq n_{2}\simeq n_{3} \simeq \frac{n}{3}$.
If the problem is easy, then $C_{1}$ has a very high accuracy, and boosting does not help.
If the problem is hard, then $C_{1}$ performs badly and $C_{2}$ will have too many samples.
Several runs of boosting are required to ensure that all of $\mathcal{D}$ is used across all component classifiers.
Each component classifier is trained with data selected according to its weight:
Samples with higher weights are more likely to be included in $D_{i+1}$.
If the samples and labels in $\mathcal{D}$ are denoted $\mathbf{x}^{i}$ and $y_{i}$, and $W_{k}(i)$ is the $k$th discrete distribution over the training, then the algorithm is…
\begin{algorithm} \caption{AdaBoost} \begin{algorithmic} \INPUT Initialize \$D\$, \$k\_\{max\}\$, \$W\_\{1\}(i)=1/n\$ for \$i=1,...,n\$, \$k = 0\$ \REPEAT \STATE \$k=k+1\$ \STATE Train \$C\_\{k\}\$ by sampling \$D\$ according to \$W\_\{k\}(i)\$ \STATE \$E\_\{k\}=\$ error of \$C\_\{k\}\$ measured on \$D\$ adjusted by \$W\_\{k\}(i)\$ \STATE \$a\_\{k\} = 0.5 *\$ ln[\$(1-E\_\{k\})/E\_\{k\}\$] \IF{\$h\_\{k\}(x^\{i\})=y\_\{i\}\$ (correct classification)} \STATE \$W\_\{k+1\} = W\_\{k\}(i) / Z\_\{k\} * e^\{-a\_\{k\}\}\$ \ELSE \STATE \$W\_\{k+1\} = W\_\{k\}(i) / Z\_\{k\} * e^\{a\_\{k\}\}\$ \ENDIF \UNTIL{\$k=k\_\{max\}\$} \RETURN \$C\_\{k\}\$ and \$a\_\{k\}\$ for \$k=1,...,k\_\{max\}\$ (i.e. the ensemble and weights) \end{algorithmic} \end{algorithm}
$E_{k}$ is the error with respect to $W_{k}(i)$; this means that $0<E_{k}<1$.
Thus $\alpha_{k}=\frac{1}{2}\ln{\left[(1-E_{k})/E_{k}\right]}$ ranges from:
If a component classifier has low error, then $\alpha_{k}$ is high (positive).
Translation: If the classifier is good, then errors are rare and we should try to add them to the training set.
If a component classifier has high error, then $\alpha_{k}$ is low (negative).
Translation: If the classifier is bad, the correct samples are more valuable, so they are weighted more.
The final classification of a point $\mathbf{x}$ is based on a weighted sum of the outputs:
$ g(\mathbf{x})=\left[\sum_{k=1}^{k_{max}}\alpha_{k}h_{k}(\mathbf{x})\right]$
The magnitude of $g(\mathbf{x})$ is the result of the different weights assigned to the component classifiers (and their signs), while the classification result is simply the sign of $g(\mathbf{x})$.
How should we set our stopping criteria? That is, how many component classifiers should we collect ($k_{max}$)?
The training error for $C_{k}$ can be written as $E_{k}=1/2-G_{k}$ for some $G_{k}>0$ (since we’re using a distribution to modulate our training error).
Then the ensemble error is simply the product:
$ E=_{k=1}^{k_{max}}\left[2\sqrt{E_{k}(1-E_{k})}\right]=_{k=1}^{k_{max}}\sqrt{1-4G_{k}^{2}}$
Thus if we keep increasing $k_{max}$, by adding more component classifiers, our error on the training set should be arbitrarily low!
“Arbitrarily low” usually means over-training.
However, this doesn’t happen often with AdaBoost.
You can’t set $k_{max}$ arbitrarily high with all of these conditions.
Nonetheless, AdaBoost is a very powerful algorithm that has been used and modified in a number of different applications.
A lot of data is partially-labeled or expensive to label.
We want to maximize label effectiveness – thus we must choose which of the unlabeled samples is most informative.
This is called “learning with queries”, “active learning”, “interactive learning”, or “cost-based learning”.
Two main types of Active Learning:
AdaBoost provides a way of combining classifiers, but there are others!
These are mixture-of-expert models, ensemble, modular, or pooled classifiers.
Assume each output is produced by a mixture model of $k$ component classifiers.
Each model has a parameter set $\boldsymbol{\theta}_{i}$.
Each estimate of the category membership for a sample $\mathbf{x}$ is $g_{ir}(\mathbf{x})=P(\omega_{r}|\mathbf{x},\boldsymbol{\theta}_{i}**$.
How do we choose $k$?
When presenting a classifier, the most basic question is:
How good is it?
There are several ways to answer this question…
Several classifiers (e.g. Bayes) gives a probability: $p(\omega_{i}|\mathbf{x})$.
We then threshold that probability, based on our risk assessment. This gives us a “hard” classification label.
Assume we have two possible classification outputs (Positive or Negative), and each of those can be right or wrong. This gives us the following table, known as a Confusion Matrix.
| Class: Positive | Class: Negative | |
|---|---|---|
| Predict: Positive | True Positive | False Positive |
| Predict: Negative | False Negative | True Negative |
From the confusion matrix, you can calculate a lot of performance metrics:
Accuracy: $\frac{TP+TN}{TP+FP+TN+FN}$
Precision (Positive Predictive Value, PPV): $\frac{TP}{TP + FP}$
Sensitivity (True Positive Rate, TPR): $\frac{TP}{TP + FN} = \frac{TP}{P}$
Specificity (True Negative Rate, TNR): $\frac{TN}{TN + FP} = \frac{TN}{N}$
Negative Predictive Value (NPV): $\frac{TN}{TN + FN}$
F1 Score: $\frac{2TP}{2TP + FP + FN}$
How can we quantify the certainty of a probabilistic classifier?
By setting our probability threshold, we can maximize sensitivity OR specificity:
Clearly we need a way to take the probability of classification into account.
The Receiver Operating Characteristics (ROC) curve is a plot of the True Positive Rate (Sensitivity) against the False Positive Rate (1 - Specificity).
The Area Under the (ROC) Curve (AUC) is a measurement of how well the classifier performs in terms of the tradeoff between the two.
Bias and variance are unavoidable; it is necessary to compensate for them.
We do this through resampling, e.g. in cross-validation.
The goal of these techniques is to determine the generalization performance.
In the simplest form, split $\mathcal{D}$ into two parts: a training and testing set.
Training set is used to create the classifier (estimate parameters, set cut-points, etc.)
After training is done, evaluate performance on the test set.
You must be careful about training on your test set! This includes:
If you have a lot of data you can use a validation set, which is used while you’re training.
This way, you can be sure that you don’t even look at your testing data until you are completely done creating your classifier.
Cross-validation is useful if your training $\mathcal{D}$ has variance, outliers, or noise.
You want to measure how your classifier changes based on which samples it trains on.
Some considerations:
If we set $k = |\mathcal{D}|$, each testing subset has one sample.
We train on $|\mathcal{D}| - 1$ samples.
What are the advantages of this approach? Disadvantages?
When the “accuracy” of a system is 90%, how sure are we that it’s not 89% or 91% or 70% instead?
Accuracy is a statistics, so it has a distribution.
We also want to calculate the mean and standard deviation of our performance metrics.
We can do this through multiple independent runs of cross-validation.
This will give a distribution of performance values, which can be reported.
Think of this as quality assurance!
In practice, $\mathcal{D}$ is far smaller than the (infinite) true dataset.
We can predict eventual accuracy using learning curves, which plot test error against the size of the training set.
These are typically described by a power-law function of the form:
$ E_{test} = a + \frac{b}{n^{\prime\alpha}} $
$ E_{train} = a - \frac{c}{n^{\prime\beta}} $
Here, $a$ is the error with a theoretically infinite sample size; ideally, it is equal to the Bayes error (optimum).
A lot of papers focus on “benchmark” datasets and evaluating new classifiers against others using a standard set of metrics.
“Is your method better than state-of-the-art (SOTA)?”
Error is not the only consideration! In some cases, reducing costs / time / resources may be just as important.
Bias, variance, and over-training must be accounted for.
Creating train / validation / testing splits is basic, standard practice for machine learning work.
We will return to a basic classification techniques: Linear Discriminants.
These are useful when your data is linearly separable, and is an introduction to more complex methods like support vector machines.
These also form the basis for neural networks, which we will discuss in the second half of the course.