Run DSE Flow¶
This page explains how to run the DSE flow and read its reports and logs.
Docker Environment (Recommend)¶
DSE infrastructure provides a Dockerfile that has all required environment ready. It is highly suggested run DSE in the docker container to enjoy the following advantages:
No worries about the DSE environment, Merlin compiler environment, and even vendor tool environment anymore.
Easy-to-use builtin commands (see instructions below).
All processes, including DSE processes, Merlin processes and vendor tool processes, are guaranteed to be terminated when the docker container is terminated.
Setup¶
Build Merlin docker image and make sure the image
melrin:latest
is available in "docker image ls".Build DSE docker image by running
./docker-build.sh
and make sure the imagemerlin-dse:latest
is available indocker image ls
.Create a new container by running
./docker/docker-run.sh /bin/bash
. Note that you need to be at the directory that contains both Merlin project and working directory when launching the container.
Execution¶
After entering to the container, you do not need to set up anything but can directly use the builtin commands as follows to make use of the DSE.
Note
You will still need to write a Python command to launch the DSE if you
would like to use more advance run-time arguments
(e.g., --disable-animation
).
See the following for detail arguments.
- Design Space Checker
This command checks if your design space definition has any problems to save your time; otherwise you may wait for several hours of DSE and found that it did nothing due to an invalid design space definition.
There have two checking modes to be used: fast and complete mode:
- fast:
Check the syntax and data type of design space definition in the config file. It usually takes a second.
- complete:
Check not only design space definition in config file but the kernel code. The errors such as missing/misspell parameter name or incorrect execution commands will be caught by this checking. Since it really runs HLS for one design point, it may take up to 30 minutes.
The docker builtin commend is used as follows:
checkds <source project dir> <working dir> <config file> <fast|complete>
- Design Space Exploration
This command launches the DSE. There have two modes to be used: fast and accurate mode:
- fast:
Perform DSE with HLS result and output a number of points with high QoR in terms of HLS cycles.
- accurate:
Perform fast mode, generates bit-streams for its outputs, and mark the one with the best QoR in terms of HLS cycle and working frequency.
The docker commend is used as follows:
autodse <source project dir> <working dir> <config file> <fast|accurate> [<database file>]
Host Machine¶
Since the DSE flow is implemented as a standard Python package, we can also execute it using Python interpreter directly:
Make sure Python 3.6+ is available.
Setup DSE package by running
python3 setup.py
. This will install dependent packages and deploy ourautodse
package in the machine. Alternatively, you can also manually install dependent packages and add the directory that containsautodse
toPYTHONPATH
.Run the design space checker to make sure the design space definition is valid:
python3 -m autodse --src-dir <source project dir> --work-dir <working dir> --config <config file> --mode <fast-check|complete-check> [--disable-animation]
Run the DSE by launching the package. Note that the optional argument
disable-animation
will use dots instead of animation to indicate the DSE is still running.python3 -m autodse --src-dir <source project dir> --work-dir <working dir> --config <config file> --db <database file> --mode <fast-dse|accurate-dse> [--disable-animation]
Outputs¶
After the DSE is done, you will find some files in the working directory:
- summary_fast.rpt
The summary file that tells you how many points have been explored in the fast mode and their results. It also provides the detail of each point so that you could know more about the search process and your design.
- summary_accurate.rpt
The summary file that tells you how many points have been processed in the accurate mode and their results. It also provides the detail of each point so that you could know more about the search process and your design.
- output
The most important directory that contains the DSE output.
- best
The symbolic link after the accurate mode pointing to the output project with the best design.
- fast
The output of fast mode. It includes a report file
output.rpt
telling you how many points we output and what are their quality. Each output is a valid Merlin project so you could continue working with it using the Merlin compiler.The output of fast mode also includes
result_dist.pdf
that draws valid HLS result on a performance-area chart. The chart also depicts a Pareto curve and shows Pareto efficient design points in the explored space.- accurate
The output of accurate mode. It includes a report file
output.rpt
telling you how many points we perform P&R and what are their quality including working frequency and actual resource utilization. The bit-stream is also available is the point could pass P&R.
- logs
All DSE log files, including the logs of overall execution flow, evaluator, explorer and the search algorithm.
- dse.log
The backup of all messages shown on the console.
- eval.log
The log of evaluation, including how many jobs have be submitted and if they encounter any issues during the evaluation.
- partX_expr.log
The log of explorer of design space partition X. It logged the time to find a better design point in the partition.
- partX_algo.log
The log of search algorithm of design space partition X. The format and contents highly depend on the search algorithm. For example,
exhaustive
simply logs all explored points;gradient
logs the process of design bottleneck analysis and important parameter identification. This log could be helpful for users to understand what DSE has done and if its exploration process makes sense or not.
- results.db
The dumped database with all evaluated results. You can simply reuse this database to continue launching the DSE if you want to explore using the fast mode again with more time.
- bak_OOOOO
The backup directory that contains all files in the working directory before launching the exploration.


Home Page
This page is part of the Angular Host app and describes the border color definitions:
-
Angular Host container
-
Angular Host router-outlet
-
Angular MFE Parent Module
-
React MFE Web Component
-
Angular MFE Static Component
-
Angular MFE Article Tile
-
Angular MFE Survey Tile


Langsdale Library
1420 Maryland Avenue, Baltimore, MD 21201
| Front Desk: |<PHONE_NUMBER> |
| Office Fax Line: |<PHONE_NUMBER> |
| Circulation: |<EMAIL_ADDRESS>| Reference: |<EMAIL_ADDRESS>|
1420 Maryland Avenue, Baltimore, MD 21201
| Front Desk: |<PHONE_NUMBER> |
| Office Fax Line: |<PHONE_NUMBER> |
| Circulation: |<EMAIL_ADDRESS>| Reference: |<EMAIL_ADDRESS>|


Name:
Hello!
{{ yourname }}
Enter your ToDo's below!
Nyron's To Do List
List
Activity Name
Add Notes
1
{{list.Errand}}
{{list.Notes}}
+


2 Kings 12:16
The trespass money and sin money was not brought into the house of the LORD: it was the priests'.


The objectives for this lesson are that by the end of the week you will be able to:
I choose this article because the media often jumps quickly on "hints" of new physics. Way too quickly in my opinion. Almost always, these hints turn out to be nothing. I am betting this particular discrepancy is probably just a fluke and I will wait until more data is collected and analyzed before getting excited.


@* * Copyright 2019 HM Revenue & Customs * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. *@ @(taxCalculationState: TaxCalculationState, previousTaxYear: Int, currentTaxYear: Int)(implicit messages: play.api.i18n.Messages) @import tags._ @heading = @{ taxCalculationState match { case _:TaxCalculationUnderpaidPaidAllState => (Messages("label.you_do_not_owe_any_more_tax"), "You do not owe any more tax") case TaxCalculationUnderpaidPaymentDueState(_, _, _, None, _) => (Messages("label.you_paid_too_little_tax_last_year"), "You paid too little tax last year") case TaxCalculationUnderpaidPaymentDueState(_, _, _, Some(_), None) => (Messages("label.you_paid_too_little_tax_last_year"), "You paid too little tax last year") case TaxCalculationUnderpaidPaymentDueState(_, _, _, Some(_), Some(SaDeadlineApproachingStatus)) => (Messages("label.you_paid_too_little_tax_last_year"), "You paid too little tax last year") case TaxCalculationUnderpaidPaymentDueState(_, _, _, Some(_), Some(SaDeadlinePassedStatus)) => (Messages("label.you_missed_the_deadline_to_pay_your_tax"), "You missed the deadline to pay your tax") case TaxCalculationUnderpaidPartPaidState(_, _, _, None, _) => (Messages("label.you_paid_too_little_tax_last_year"), "You paid too little tax last year") case TaxCalculationUnderpaidPartPaidState(_, _, _, Some(_), None) => (Messages("label.you_paid_too_little_tax_last_year"), "You paid too little tax last year") case TaxCalculationUnderpaidPartPaidState(_, _, _, Some(_), Some(SaDeadlineApproachingStatus)) => (Messages("label.you_paid_too_little_tax_last_year"), "You paid too little tax last year") case TaxCalculationUnderpaidPartPaidState(_, _, _, Some(_), Some(SaDeadlinePassedStatus)) => (Messages("label.you_missed_the_deadline_to_pay_your_tax"), "You missed the deadline to pay your tax") case _:TaxCalculationOverpaidState => (Messages("label.you_paid_too_much_tax_last_year"), "You paid too much tax last year") case _:TaxCalculationDisabledState => (Messages("label.your_income_tax_for_last_year"), "Your Income Tax for last year") case _ => ("", "") } } @contentAndLinks = @{ taxCalculationState match { case TaxCalculationOverpaidRefundState(amount, startOfTaxYear, endOfTaxYear) => (Messages("label.hmrc_owes_you_a_refund_for_the_to_tax_year", amount, startOfTaxYear.toString, endOfTaxYear.toString), List( (Messages("label.claim_your_tax_refund"), "/tax-you-paid/status", "Claim your tax refund"), (Messages("label.find_out_why_you_paid_too_much"), s"/tax-you-paid/$previousTaxYear-$currentTaxYear/paid-too-much/reasons", "Find out why you paid too much") )) case TaxCalculationOverpaidPaymentProcessingState(amount) => (Messages("label.hmrc_is_processing_your_refund", amount), List( (Messages("label.find_out_why_you_paid_too_much"), s"/tax-you-paid/$previousTaxYear-$currentTaxYear/paid-too-much/reasons", "Find out why you paid too much") )) case TaxCalculationOverpaidPaymentPaidState(amount, datePaid) => (Messages("label.hmrc_has_paid_your_refund", amount), List( (Messages("label.find_out_why_you_paid_too_much"), s"/tax-you-paid/$previousTaxYear-$currentTaxYear/paid-too-much/reasons", "Find out why you paid too much") )) case TaxCalculationOverpaidPaymentChequeSentState(amount, datePaid) => (Messages("label.hmrc_sent_you_a_cheque_for", amount), List( (Messages("label.find_out_why_you_paid_too_much"), s"/tax-you-paid/$previousTaxYear-$currentTaxYear/paid-too-much/reasons", "Find out why you paid too much") )) case TaxCalculationUnderpaidPaymentDueState(amount, startOfTaxYear, endOfTaxYear, None, None) => (Messages("label.you_owe_hmrc_for_the_to_tax_year", amount, startOfTaxYear.toString, endOfTaxYear.toString), List( (Messages("label.find_out_why_you_paid_too_little"), s"/tax-you-paid/$previousTaxYear-$currentTaxYear/paid-too-little/reasons", "Find out why you paid too little") )) case TaxCalculationUnderpaidPaymentDueState(amount, startOfTaxYear, endOfTaxYear, Some(dueDate), None) => (Messages("label.you_owe_hmrc_for_the_to_tax_year_you_must_pay_by_", amount, startOfTaxYear.toString, endOfTaxYear.toString, dueDate), List( (Messages("label.make_a_payment"), "https://www.gov.uk/simple-assessment", "Make a payment"), (Messages("label.find_out_why_you_paid_too_little"), s"/tax-you-paid/$previousTaxYear-$currentTaxYear/paid-too-little/reasons", "Find out why you paid too little") )) case TaxCalculationUnderpaidPaymentDueState(amount, startOfTaxYear, endOfTaxYear, Some(dueDate), Some(saDeadlineStatus)) => saDeadlineStatus match { case SaDeadlineApproachingStatus => (Messages("label.you_owe_hmrc_for_the_to_tax_year_you_must_pay_by_", amount, startOfTaxYear.toString, endOfTaxYear.toString, dueDate), List( (Messages("label.make_a_payment"), "https://www.gov.uk/simple-assessment", "Make a payment"), (Messages("label.find_out_why_you_paid_too_little"), s"/tax-you-paid/$previousTaxYear-$currentTaxYear/paid-too-little/reasons", "Find out why you paid too little") )) case SaDeadlinePassedStatus => (Messages("label.you_owe_hmrc_for_the_to_tax_year_you_should_have_paid_", amount, startOfTaxYear.toString, endOfTaxYear.toString, dueDate), List( (Messages("label.make_a_payment"), "https://www.gov.uk/simple-assessment", "You missed the deadline to pay your tax") )) } case TaxCalculationUnderpaidPartPaidState(amount, startOfTaxYear, endOfTaxYear, None, None) => (Messages("label.you_owe_hmrc_for_the_to_tax_year", amount, startOfTaxYear.toString, endOfTaxYear.toString), List( (Messages("label.find_out_why_you_paid_too_little"), s"/tax-you-paid/$previousTaxYear-$currentTaxYear/paid-too-little/reasons", "Find out why you paid too little") )) case TaxCalculationUnderpaidPartPaidState(amount, startOfTaxYear, endOfTaxYear, Some(dueDate), None) => (Messages("label.you_still_owe_hmrc_for_the_to_tax_year_you_must_pay_by_", amount, startOfTaxYear.toString, endOfTaxYear.toString, dueDate), List( (Messages("label.make_a_payment"), "https://www.gov.uk/simple-assessment", "Make a payment"), (Messages("label.find_out_why_you_paid_too_little"), s"/tax-you-paid/$previousTaxYear-$currentTaxYear/paid-too-little/reasons", "Find out why you paid too little") )) case TaxCalculationUnderpaidPartPaidState(amount, startOfTaxYear, endOfTaxYear, Some(dueDate), Some(saDeadlineStatus)) => saDeadlineStatus match { case SaDeadlineApproachingStatus => (Messages("label.you_still_owe_hmrc_for_the_to_tax_year_you_must_pay_by_", amount, startOfTaxYear.toString, endOfTaxYear.toString, dueDate), List( (Messages("label.make_a_payment"), "https://www.gov.uk/simple-assessment", "Make a payment"), (Messages("label.find_out_why_you_paid_too_little"), s"/tax-you-paid/$previousTaxYear-$currentTaxYear/paid-too-little/reasons", "Find out why you paid too little") )) case SaDeadlinePassedStatus => { (Messages("label.you_still_owe_hmrc_for_the_to_tax_year_you_should_have_paid_", amount, startOfTaxYear.toString, endOfTaxYear.toString, dueDate), List( (Messages("label.make_a_payment"), "https://www.gov.uk/simple-assessment", "You missed the deadline to pay your tax") )) } } case TaxCalculationUnderpaidPaidAllState(startOfTaxYear, endOfTaxYear, _) => (Messages("label.you_have_no_payments_to_make_to_hmrc_for_the_to_tax_year", startOfTaxYear.toString, endOfTaxYear.toString), List( (Messages("label.find_out_why_you_paid_too_little"), s"/tax-you-paid/$previousTaxYear-$currentTaxYear/paid-too-little/reasons", "Find out why you paid too little") )) case TaxCalculationDisabledState(startOfTaxYear, endOfTaxYear) => (Messages("label.check_to_see_if_you_paid_the_right_amount_of_tax_from_", startOfTaxYear.toString, endOfTaxYear.toString), List()) case _ => ("", List()) } } @card( url = Some("/tax-you-paid/status"), gaAction = Some("Income"), gaLabel = Some(heading._2), heading = heading._1, headingTag = "h3", bodyContent = Some(Html("
" + contentAndLinks._1 + "
")) ) { @if(!contentAndLinks._2.isEmpty){


15 GEE, Continuous Outcome: Beat the Blues
15.1 Packages
15.1.1 CRAN
library(tidyverse) # all things tidy
library(pander) # nice looking genderal tabulations
library(furniture) # nice table1() descriptives
library(texreg) # Convert Regression Output to LaTeX or HTML Tables
library(psych) # contains some useful functions, like headTail
library(interactions)
library(performance)
library(lme4) # Linear, generalized linear, & nonlinear mixed models
library(corrplot) # Vizualize correlation matrix
library(gee) # Genderalized Estimation Equation Solver
library(geepack) # Genderalized Estimation Equation Package
library(MuMIn) # Multi-Model Inference (caluclate QIC)
library(HSAUR) # package with the dataset
15.2 Background
This dataset was used as an example in Chapter 11 of “A Handbook of Statistical Analysis using R” by Brian S. Everitt and Torsten Hothorn. The authors include this data set in their HSAUR package on
CRAN
.
The Background
“These data resulted from a clinical trial of an interactive multimedia program called ‘Beat the Blues’ designed to deliver cognitive behavioral therapy (CBT) to depressed patients via a computer terminal. Full details are given here: Proudfoot et. al (2003), but in essence Beat the Blues is an interactive program using multimedia techniques, in particular video vignettes.”
“In a randomized controlled trial of the program, patients with depression recruited in primary care were randomized to either the Beating the Blues program or to”Treatment as Usual” (TAU). Patients randomized to the BEat the Blues also received pharmacology and/or general practice (GP) support and practical/social help,offered as part of treatment as usual, with the exception of any face-to-face counseling or psychological intervention. Patients allocated to TAU received whatever treatment their GP prescribed. The latter included, besides any medication, discussion of problems with GP, provisions of practical/social help, referral to a counselor, referral to a practice nurse, referral to mental health professionals, or further physical examination.
The Research Question
Net of other factors (use of antidepressants and length of the current episode), does the Beat-the-Blues program results in better depression trajectories over treatment as usual?
The Data
The variables are as follows:
drug
did the patient take anti-depressant drugs (No or Yes)
length
the length of the current episode of depression, a factor with levels:- “<6m” less than six months
- “>6m” more than six months
- “<6m” less than six months
treatment
treatment group, a factor with levels:- “TAU” treatment as usual
- “BtheB” Beat the Blues
- “TAU” treatment as usual
bdi.pre
Beck Depression Inventory II, before treatment
bdi.2m
Beck Depression Inventory II, after 2 months
bdi.4m
Beck Depression Inventory II, after 4 months
bdi.6m
Beck Depression Inventory II, after 6 months
bdi.8m
Beck Depression Inventory II, after 8 months
15.2.1 Read in the data
data(BtheB, package = "HSAUR")
%>%
BtheB ::headTail() psych
drug length treatment bdi.pre bdi.2m bdi.4m bdi.6m bdi.8m
1 No >6m TAU 29 2 2 <NA> <NA>
2 Yes >6m BtheB 32 16 24 17 20
3 Yes <6m TAU 25 20 <NA> <NA> <NA>
4 No >6m BtheB 21 17 16 10 9
... <NA> <NA> <NA> ... ... ... ... ...
97 Yes <6m TAU 28 <NA> <NA> <NA> <NA>
98 No >6m BtheB 11 22 9 11 11
99 No <6m TAU 13 5 5 0 6
100 Yes <6m TAU 43 <NA> <NA> <NA> <NA>
15.2.2 Wide Format
Tidy up the dataset
WIDE format * One line per person * Good for descriptive tables & correlation matrices
<- BtheB %>%
btb_wide ::mutate(id = row_number()) %>% # create a new variable to ID participants
dplyr::select(id, treatment, # specify that ID variable is first
dplyr
drug, length,.2m, bdi.4m, bdi.6m, bdi.8m) bdi.pre, bdi
%>%
btb_wide ::headTail(top = 10, bottom = 10) psych
id treatment drug length bdi.pre bdi.2m bdi.4m bdi.6m bdi.8m
1 1 TAU No >6m 29 2 2 <NA> <NA>
2 2 BtheB Yes >6m 32 16 24 17 20
3 3 TAU Yes <6m 25 20 <NA> <NA> <NA>
4 4 BtheB No >6m 21 17 16 10 9
5 5 BtheB Yes >6m 26 23 <NA> <NA> <NA>
6 6 BtheB Yes <6m 7 0 0 0 0
7 7 TAU Yes <6m 17 7 7 3 7
8 8 TAU No >6m 20 20 21 19 13
9 9 BtheB Yes <6m 18 13 14 20 11
10 10 BtheB Yes >6m 20 5 5 8 12
... ... <NA> <NA> <NA> ... ... ... ... ...
91 91 TAU No <6m 16 <NA> <NA> <NA> <NA>
92 92 BtheB Yes <6m 30 26 28 <NA> <NA>
93 93 BtheB Yes <6m 17 8 7 12 <NA>
94 94 BtheB No >6m 19 4 3 3 3
95 95 BtheB No >6m 16 11 4 2 3
96 96 BtheB Yes >6m 16 16 10 10 8
97 97 TAU Yes <6m 28 <NA> <NA> <NA> <NA>
98 98 BtheB No >6m 11 22 9 11 11
99 99 TAU No <6m 13 5 5 0 6
100 100 TAU Yes <6m 43 <NA> <NA> <NA> <NA>
15.2.3 Long Format
Restructure to long format
LONG FORMAT * One line per observation * Good for plots and modeling
<- btb_wide %>%
btb_long ::pivot_longer(cols = c(bdi.2m, bdi.4m, bdi.6m, bdi.8m), # all existing variables (not quoted)
tidyrnames_to = "month",
names_pattern = "bdi.(.)m",
values_to = "bdi") %>%
::mutate(month = as.numeric(month)) %>%
dplyr::filter(complete.cases(id, bdi, treatment, month)) %>%
dplyr::arrange(id, month) %>%
dplyr::select(id, treatment, drug, length, bdi.pre, month, bdi) dplyr
%>%
btb_long ::headTail(top = 10, bottom = 10) psych
id treatment drug length bdi.pre month bdi
1 1 TAU No >6m 29 2 2
2 1 TAU No >6m 29 4 2
3 2 BtheB Yes >6m 32 2 16
4 2 BtheB Yes >6m 32 4 24
5 2 BtheB Yes >6m 32 6 17
6 2 BtheB Yes >6m 32 8 20
7 3 TAU Yes <6m 25 2 20
8 4 BtheB No >6m 21 2 17
9 4 BtheB No >6m 21 4 16
10 4 BtheB No >6m 21 6 10
11 ... <NA> <NA> <NA> ... ... ...
12 96 BtheB Yes >6m 16 6 10
13 96 BtheB Yes >6m 16 8 8
14 98 BtheB No >6m 11 2 22
15 98 BtheB No >6m 11 4 9
16 98 BtheB No >6m 11 6 11
17 98 BtheB No >6m 11 8 11
18 99 TAU No <6m 13 2 5
19 99 TAU No <6m 13 4 5
20 99 TAU No <6m 13 6 0
21 99 TAU No <6m 13 8 6
15.3 Exploratory Data Analysis
15.3.1 Visualize: Person-profile Plots
Create spaghetti plots of the raw, observed data
%>%
btb_long ggplot(aes(x = month,
y = bdi)) +
geom_point() +
geom_line(aes(group = id),
size = 1,
alpha = 0.3) +
geom_smooth(method = "lm") +
theme_bw() +
facet_grid(.~ treatment) +
labs(title = "BtheB - Observed Data Across Time with LM Smoother",
subtitle = "Seperate by Treatment")
%>%
btb_long ggplot(aes(x = month,
y = bdi)) +
geom_point() +
geom_line(aes(group = id),
size = 1,
alpha = 0.3) +
geom_smooth(method = "lm") +
facet_grid(drug~ treatment, labeller = label_both) +
theme_bw() +
labs(title = "BtheB - Observed Data Across Time with LM Smoother",
subtitle = "Seperate by Treatment & Antidepressant Use")
%>%
btb_long ggplot(aes(x = month,
y = bdi)) +
geom_point() +
geom_line(aes(group = id,
color = length),
size = 1,
alpha = 0.3) +
geom_smooth(aes(color = length),
method = "lm",
size = 1.25,
se = FALSE) +
facet_grid(drug~ treatment, labeller = label_both) +
theme_bw() +
labs(title = "BtheB - Observed Data Across Time with LM Smoother",
subtitle = "Seperate by Treatment, Antidepressant Use, & Length of Current Episode")
%>%
btb_long ggplot(aes(x = month,
y = bdi,
color = treatment,
fill = treatment)) +
geom_smooth(method = "lm") +
theme_bw() +
labs(title = "BtheB - Predictions from TWO Seperate Single Simple Linear Models (lm)",
subtitle = "Assumes Independence of Repeated Measures")
15.3.2 Calculate the Observed Correlation Structure
<- btb_wide %>%
bdi_corr ::select(starts_with("bdi")) %>%
dplyr::cor(use="pairwise.complete.obs")
stats
bdi_corr
bdi.pre bdi.2m bdi.4m bdi.6m bdi.8m
bdi.pre 1.0000000 0.6142207 0.5691248 0.5077286 0.3835090
bdi.2m 0.6142207 1.0000000 0.7903346 0.7849188 0.7038158
bdi.4m 0.5691248 0.7903346 1.0000000 0.8166591 0.7220149
bdi.6m 0.5077286 0.7849188 0.8166591 1.0000000 0.8107773
bdi.8m 0.3835090 0.7038158 0.7220149 0.8107773 1.0000000
15.4 Multiple Regression (OLS)
This ignores any correlation between repeated measures on the same individual and treats all observations as independent.
15.4.1 Fit the models
<- stats::lm(bdi ~ bdi.pre + length + drug + treatment + month,
btb_lm_1 data = btb_long)
<- stats::lm(bdi ~ bdi.pre + length + drug + treatment*month,
btb_lm_2 data = btb_long)
<- stats::lm(bdi ~ bdi.pre + length + drug + treatment + drug*month,
btb_lm_3 data = btb_long)
<- stats::lm(bdi ~ bdi.pre + length + drug*treatment*month,
btb_lm_4 data = btb_long)
15.4.2 Parameter Estimates Table
::knitreg(list(btb_lm_1, btb_lm_2, btb_lm_3, btb_lm_4),
texreglabel = "lm",
caption = "OLS")
| Model 1 | Model 2 | Model 3 | Model 4 | |
|---|---|---|---|---|
| (Intercept) | 7.88*** | 7.77*** | 7.21*** | 7.33** |
| (1.78) | (2.08) | (2.03) | (2.30) | |
| bdi.pre | 0.57*** | 0.57*** | 0.57*** | 0.56*** |
| (0.05) | (0.05) | (0.05) | (0.05) | |
| length>6m | 1.75 | 1.75 | 1.78 | 1.86 |
| (1.11) | (1.11) | (1.11) | (1.10) | |
| drugYes | -3.55** | -3.55** | -2.10 | -2.00 |
| (1.14) | (1.15) | (2.39) | (3.75) | |
| treatmentBtheB | -3.35** | -3.13 | -3.36** | -3.31 |
| (1.10) | (2.36) | (1.10) | (3.13) | |
| month | -0.96*** | -0.93** | -0.82** | -0.60 |
| (0.23) | (0.34) | (0.31) | (0.40) | |
| treatmentBtheB:month | -0.05 | -0.56 | ||
| (0.47) | (0.63) | |||
| drugYes:month | -0.32 | -1.02 | ||
| (0.47) | (0.73) | |||
| drugYes:treatmentBtheB | -0.23 | |||
| (4.92) | ||||
| drugYes:treatmentBtheB:month | 1.31 | |||
| (0.98) | ||||
| R2 | 0.40 | 0.40 | 0.40 | 0.42 |
| Adj. R2 | 0.39 | 0.38 | 0.39 | 0.40 |
| Num. obs. | 280 | 280 | 280 | 280 |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
15.4.3 Plot the model predictions
::interact_plot(model = btb_lm_1,
interactionspred = month,
modx = treatment,
main.title = "LM with 95% CI",
interval = TRUE) # default = 95% confidence interval
::interact_plot(model = btb_lm_1,
interactionspred = month,
modx = treatment,
main.title = "LM with 68% CI (M +/- SEM)",
interval = TRUE,
int.width = 0.68) # change to 68% CI, which is +/- 1 SEM
::interact_plot(model = btb_lm_1,
interactionspred = month,
modx = treatment,
interval = TRUE,
colors = rep("black", times = 3),
x.label = "Month",
y.label = "Predicted BDI",
legend.main = "Baseline BDI:") +
theme_bw() +
theme(legend.position = "bottom",
legend.key.width = unit(2, "cm")) +
labs(title = "BtheB - Predictions from the LM model (95% CI)",
subtitle = "Trajectory for a person with BL depression < 6 months and randomized to TAU")
::Effect(focal.predictors = c("treatment", "month"),
effectsmod = btb_lm_1) %>%
%>%
data.frame ::mutate(treatment = fct_reorder2(treatment, month, fit)) %>%
dplyrggplot(aes(x = month,
y = fit)) +
geom_line(aes(color = treatment)) +
geom_ribbon(aes(ymin = lower, # 95% confidence intervals
ymax = upper,
fill = treatment),
alpha = 0.3) +
geom_ribbon(aes(ymin = fit - se, # Mean +/- 1 standard error for the mean (SEM)
ymax = fit + se,
fill = treatment),
alpha = 0.3) +
theme_bw() +
labs(title = "BtheB - Predictions from a Single Linear Model (lm)",
subtitle = "Assumes Independence of Repeated Measures") +
theme(legend.position = c(1, 1),
legend.justification = c(1.1, 1.1),
legend.background = element_rect(color = "black"))
15.5 Multilevel Models (MLM)
15.5.1 Fit the models
<- lmerTest::lmer(bdi ~ bdi.pre + length + drug + treatment + month + (1 | id),
btb_lmer_RI data = btb_long,
REML = TRUE)
<- lmerTest::lmer(bdi ~ bdi.pre + length + drug + treatment + month + (month | id),
btb_lmer_RIAS data = btb_long,
REML = TRUE,
control = lmerControl(optimizer = "Nelder_Mead"))
15.5.2 Parameter Estimates Table
::knitreg(list(btb_lm_1, btb_lmer_RI, btb_lmer_RIAS),
texregcustom.model.names = c("OLS", "MLM-RI", "MLM-RIAS"),
label = "mlm",
caption = "LM vs. MLM")
| OLS | MLM-RI | MLM-RIAS | |
|---|---|---|---|
| (Intercept) | 7.88*** | 5.92* | 5.94** |
| (1.78) | (2.31) | (2.30) | |
| bdi.pre | 0.57*** | 0.64*** | 0.64*** |
| (0.05) | (0.08) | (0.08) | |
| length>6m | 1.75 | 0.24 | 0.10 |
| (1.11) | (1.68) | (1.67) | |
| drugYes | -3.55** | -2.79 | -2.89 |
| (1.14) | (1.77) | (1.76) | |
| treatmentBtheB | -3.35** | -2.36 | -2.49 |
| (1.10) | (1.71) | (1.71) | |
| month | -0.96*** | -0.71*** | -0.70*** |
| (0.23) | (0.15) | (0.16) | |
| R2 | 0.40 | ||
| Adj. R2 | 0.39 | ||
| Num. obs. | 280 | 280 | 280 |
| AIC | 1882.08 | 1885.16 | |
| BIC | 1911.16 | 1921.50 | |
| Log Likelihood | -933.04 | -932.58 | |
| Num. groups: id | 97 | 97 | |
| Var: id (Intercept) | 51.44 | 50.56 | |
| Var: Residual | 25.27 | 23.87 | |
| Var: id month | 0.23 | ||
| Cov: id (Intercept) month | -0.31 | ||
| ***p < 0.001; **p < 0.01; *p < 0.05 |
15.5.3 Likelihood Ratio Test
anova(btb_lmer_RI,
btb_lmer_RIAS, refit = FALSE)
Data: btb_long
Models:
btb_lmer_RI: bdi ~ bdi.pre + length + drug + treatment + month + (1 | id)
btb_lmer_RIAS: bdi ~ bdi.pre + length + drug + treatment + month + (month | id)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
btb_lmer_RI 8 1882.1 1911.2 -933.04 1866.1
btb_lmer_RIAS 10 1885.2 1921.5 -932.58 1865.2 0.9236 2 0.6301
::compare_performance(btb_lm_1,
performance
btb_lmer_RI,
btb_lmer_RIAS,rank = TRUE)
# Comparison of Model Performance Indices
Name | Model | RMSE | Sigma | AIC weights | BIC weights | Performance-Score
-----------------------------------------------------------------------------------------------
btb_lmer_RI | lmerModLmerTest | 4.234 | 5.027 | 0.831 | 0.995 | 97.86%
btb_lmer_RIAS | lmerModLmerTest | 4.016 | 4.885 | 0.169 | 0.005 | 55.20%
btb_lm_1 | lm | 8.560 | 8.654 | 8.42e-28 | 6.20e-27 | 0.00%
15.5.4 Plot the model predictions
::interact_plot(model = btb_lmer_RI,
interactionspred = month,
modx = treatment,
mod2 = drug,
interval = TRUE) +
theme_bw() +
labs(title = "BtheB - Predictions from a Multilevel Model (lmer)") +
theme(legend.position = c(0, 0),
legend.justification = c(-0.1, -0.1),
legend.background = element_rect(color = "black"))
::Effect(c("treatment", "month", "drug"),
effectsmod = btb_lmer_RI) %>%
%>%
data.frame ::mutate(treatment = fct_reorder2(treatment, month, fit)) %>%
dplyrggplot(aes(x = month,
y = fit)) +
geom_line(aes(color = treatment)) +
geom_ribbon(aes(ymin = lower,
ymax = upper,
fill = treatment),
alpha = 0.3) +
theme_bw() +
facet_grid(.~ drug, labeller = label_both) +
labs(title = "BtheB - Predictions from a Multilevel Model (lmer)") +
theme(legend.position = c(0, 0),
legend.justification = c(-0.1, -0.1),
legend.background = element_rect(color = "black"))
15.6 General Estimating Equations, GEE
15.6.1 Explore Various Correlation Structures
<IP_ADDRESS> Fit the models - Main effects to determine correlation structure
Use the gee()
function from the gee
package for the results to be used in a texreg::knitreg()
tables.
The output below each model is the ‘starting’ model assuming independence, so they will all be the same here.
<- gee::gee(bdi ~ bdi.pre + length + drug + treatment + month,
btb_gee_in data = btb_long,
id = id,
family = gaussian,
corstr = 'independence')
(Intercept) bdi.pre length>6m drugYes treatmentBtheB
7.8830747 0.5723729 1.7530800 -3.5460058 -3.3539662
month
-0.9608077
<- gee::gee(bdi ~ bdi.pre + length + drug + treatment + month,
btb_gee_ex data = btb_long,
id = id,
family = gaussian,
corstr = 'exchangeable')
(Intercept) bdi.pre length>6m drugYes treatmentBtheB
7.8830747 0.5723729 1.7530800 -3.5460058 -3.3539662
month
-0.9608077
# The AR-1 fails if any subjects have only 1 observation
# to use this one, we would need to remove participants with only 1 BDI
# btb_gee_ar <- gee(bdi ~ bdi.pre + length + drug + treatment + month,
# data = btb_long,
# id = id,
# family = gaussian,
# corstr = 'AR-M',
# Mv = 1)
<- gee::gee(bdi ~ bdi.pre + length + drug + treatment + month,
btb_gee_un data = btb_long,
id = id,
family = gaussian,
corstr = 'unstructured')
(Intercept) bdi.pre length>6m drugYes treatmentBtheB
7.8830747 0.5723729 1.7530800 -3.5460058 -3.3539662
month
-0.9608077
summary(btb_gee_in)
GEE: GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
gee S-function, version 4.13 modified 98/01/27 (1998)
Model:
Link: Identity
Variance to Mean Relation: Gaussian
Correlation Structure: Independent
Call:
gee::gee(formula = bdi ~ bdi.pre + length + drug + treatment +
month, id = id, data = btb_long, family = gaussian, corstr = "independence")
Summary of Residuals:
Min 1Q Median 3Q Max
-24.20158432 -5.31202378 0.01101526 5.29503741 27.77789553
Coefficients:
Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept) 7.8830747 1.78048908 4.427477 2.19973944 3.583640
bdi.pre 0.5723729 0.05486079 10.433188 0.08853253 6.465114
length>6m 1.7530800 1.10849861 1.581490 1.41954159 1.234962
drugYes -3.5460058 1.14469086 -3.097785 1.73069664 -2.048889
treatmentBtheB -3.3539662 1.09831939 -3.053726 1.71390982 -1.956909
month -0.9608077 0.23263437 -4.130119 0.17688635 -5.431780
Estimated Scale Parameter: 74.8854
Number of Iterations: 1
Working Correlation
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 0 1 0
[4,] 0 0 0 1
summary(btb_gee_ex)
GEE: GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
gee S-function, version 4.13 modified 98/01/27 (1998)
Model:
Link: Identity
Variance to Mean Relation: Gaussian
Correlation Structure: Exchangeable
Call:
gee::gee(formula = bdi ~ bdi.pre + length + drug + treatment +
month, id = id, data = btb_long, family = gaussian, corstr = "exchangeable")
Summary of Residuals:
Min 1Q Median 3Q Max
-25.4478843 -6.3276726 -0.8152833 4.3622258 25.4078115
Coefficients:
Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept) 5.8855129 2.32380381 2.5327065 2.10712166 2.7931529
bdi.pre 0.6399964 0.08033495 7.9665999 0.07931263 8.0692874
length>6m 0.2084783 1.69179766 0.1232288 1.48052530 0.1408137
drugYes -2.7742506 1.78397557 -1.5550945 1.64824318 -1.6831561
treatmentBtheB -2.3360241 1.72621751 -1.3532617 1.66217026 -1.4054060
month -0.7078407 0.14254124 -4.9658660 0.15394156 -4.5981134
Estimated Scale Parameter: 77.14393
Number of Iterations: 5
Working Correlation
[,1] [,2] [,3] [,4]
[1,] 1.0000000 0.6915241 0.6915241 0.6915241
[2,] 0.6915241 1.0000000 0.6915241 0.6915241
[3,] 0.6915241 0.6915241 1.0000000 0.6915241
[4,] 0.6915241 0.6915241 0.6915241 1.0000000
summary(btb_gee_un)
GEE: GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
gee S-function, version 4.13 modified 98/01/27 (1998)
Model:
Link: Identity
Variance to Mean Relation: Gaussian
Correlation Structure: Unstructured
Call:
gee::gee(formula = bdi ~ bdi.pre + length + drug + treatment +
month, id = id, data = btb_long, family = gaussian, corstr = "unstructured")
Summary of Residuals:
Min 1Q Median 3Q Max
-25.1527937 -6.1091139 -0.5896205 4.7316139 25.9041542
Coefficients:
Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept) 6.3905215 2.28769760 2.793429 2.15668950 2.9631162
bdi.pre 0.6171798 0.07744569 7.969195 0.08081777 7.6366846
length>6m 0.5834398 1.61626647 0.360980 1.46837275 0.3973377
drugYes -2.7908835 1.69816226 -1.643473 1.63741987 -1.7044398
treatmentBtheB -2.4261698 1.64272613 -1.476917 1.65519523 -1.4657907
month -0.7628336 0.18121518 -4.209546 0.15643591 -4.8763329
Estimated Scale Parameter: 76.40371
Number of Iterations: 5
Working Correlation
[,1] [,2] [,3] [,4]
[1,] 1.0000000 0.7069560 0.5704892 0.4714744
[2,] 0.7069560 1.0000000 0.6086188 0.4637445
[3,] 0.5704892 0.6086188 1.0000000 0.5454963
[4,] 0.4714744 0.4637445 0.5454963 1.0000000
<IP_ADDRESS> Parameter Estimates Table
::knitreg(list(btb_lm_1,
texreg
btb_lmer_RI,
btb_gee_in,
btb_gee_ex,
btb_gee_un),custom.model.names = c("OLS",
"MLM-RI",
"GEE-in",
"GEE-ex",
"GEE-un"),
label = "GEEs",
caption = "LM, MLM, and GEE")
| OLS | MLM-RI | GEE-in | GEE-ex | GEE-un | |
|---|---|---|---|---|---|
| (Intercept) | 7.88*** | 5.92* | 7.88*** | 5.89** | 6.39** |
| (1.78) | (2.31) | (2.20) | (2.11) | (2.16) | |
| bdi.pre | 0.57*** | 0.64*** | 0.57*** | 0.64*** | 0.62*** |
| (0.05) | (0.08) | (0.09) | (0.08) | (0.08) | |
| length>6m | 1.75 | 0.24 | 1.75 | 0.21 | 0.58 |
| (1.11) | (1.68) | (1.42) | (1.48) | (1.47) | |
| drugYes | -3.55** | -2.79 | -3.55* | -2.77 | -2.79 |
| (1.14) | (1.77) | (1.73) | (1.65) | (1.64) | |
| treatmentBtheB | -3.35** | -2.36 | -3.35 | -2.34 | -2.43 |
| (1.10) | (1.71) | (1.71) | (1.66) | (1.66) | |
| month | -0.96*** | -0.71*** | -0.96*** | -0.71*** | -0.76*** |
| (0.23) | (0.15) | (0.18) | (0.15) | (0.16) | |
| R2 | 0.40 | ||||
| Adj. R2 | 0.39 | ||||
| Num. obs. | 280 | 280 | 280 | 280 | 280 |
| AIC | 1882.08 | ||||
| BIC | 1911.16 | ||||
| Log Likelihood | -933.04 | ||||
| Num. groups: id | 97 | ||||
| Var: id (Intercept) | 51.44 | ||||
| Var: Residual | 25.27 | ||||
| Scale | 74.89 | 77.14 | 76.40 | ||
| ***p < 0.001; **p < 0.01; *p < 0.05 |
<IP_ADDRESS> Re-Fit Models
Use the geeglm()
function from the geepack
package for the results to be used in a anova()
table and interaction plots.
This function does NOT produce the same starting model output as gee::gee()
.
<- geepack::geeglm(bdi ~ bdi.pre + length + drug + treatment + month,
btb_geeglm_in data = btb_long,
id = id,
wave = month,
family = gaussian,
corstr = 'independence')
<- geepack::geeglm(bdi ~ bdi.pre + length + drug + treatment + month,
btb_geeglm_ex data = btb_long,
id = id,
wave = month,
family = gaussian,
corstr = 'exchangeable')
<- geepack::geeglm(bdi ~ bdi.pre + length + drug + treatment + month,
btb_geeglm_ar data = btb_long,
id = id,
wave = month,
family = gaussian,
corstr = 'ar1')
<- geepack::geeglm(bdi ~ bdi.pre + length + drug + treatment + month,
btb_geeglm_un data = btb_long,
id = id,
wave = month,
family = gaussian,
corstr = 'unstructured')
Can’t Use the Likelihood Ratio Test
The anova()
function is used to compare nested models for parameters (fixed effects), not correlation structures.
anova(btb_geeglm_in, btb_geeglm_ex)
Models are identical
NULL
anova(btb_geeglm_in, btb_geeglm_ar)
Models are identical
NULL
anova(btb_geeglm_in, btb_geeglm_un)
Models are identical
NULL
<IP_ADDRESS> Variaous QIC Measures of Fit
References:
Pan, W. 2001. Akaike’s information criterion in generalized estimating equations. Biometrics 57:120-125. https://onlinelibrary.wiley.com/doi/abs/10.1111/j.0006-341X.2001.00120.x
Burnham, K. P. and D. R. Anderson. 2002. Model selection and multimodel inference: a practical information-theoretic approach. Second edition. Springer Science and Business Media, Inc., New York. https://cds.cern.ch/record/1608735/files/9780387953649_TOC.pdf
The QIC()
is one way to try to measure model fit. You can enter more than one model into a single function call.
QIC(I) based on independence model <– suggested by Pan (Biometric, March 2001), asymptotically unbiased estimator (choose the correlation stucture that produces the smallest QIC(I), p122)
::QIC(btb_geeglm_in,
MuMIn
btb_geeglm_ex,
btb_geeglm_ar,
btb_geeglm_un, typeR = FALSE) %>% # default
::pander(caption = "QIC") pander
| QIC | |
|---|---|
| btb_geeglm_in | 307 |
| btb_geeglm_ex | 296 |
| btb_geeglm_ar | 298 |
| btb_geeglm_un | 297 |
QIC(R) is based on quasi-likelihood of a working correlation R model, can NOT be used to select the working correlation matrix.
::QIC(btb_geeglm_in,
MuMIn
btb_geeglm_ex,
btb_geeglm_ar,
btb_geeglm_un, typeR = TRUE) # NOT the default
QIC
btb_geeglm_in 306.5589
btb_geeglm_ex 304.5003
btb_geeglm_ar 304.6425
btb_geeglm_un 304.4087
QIC_U(R) approximates QIC(R), and while both are useful for variable selection, they can NOT be applied to select the working correlation matrix.
::QICu(btb_geeglm_in,
MuMIn
btb_geeglm_ex,
btb_geeglm_ar, btb_geeglm_un)
QICu
btb_geeglm_in 292.0000
btb_geeglm_ex 283.7551
btb_geeglm_ar 285.6132
btb_geeglm_un 284.1707
::model.sel(btb_geeglm_in,
MuMIn
btb_geeglm_ex,
btb_geeglm_ar,
btb_geeglm_un, rank = "QIC") #sorts the best to the TOP, uses QIC(I)
Model selection table
(Int) bdi.pre drg lng mnt trt corstr qLik QIC delta weight
btb_geeglm_ex 5.880 0.6402 + + -0.7070 + exchng -140 296.3 0.00 0.450
btb_geeglm_un 6.068 0.6307 + + -0.7061 + unstrc -140 296.6 0.32 0.382
btb_geeglm_ar 6.620 0.5956 + + -0.7357 + ar1 -140 298.3 2.00 0.165
btb_geeglm_in 7.883 0.5724 + + -0.9608 + indpnd -140 306.6 10.30 0.003
Abbreviations:
corstr: exchng = 'exchangeable', indpnd = 'independence',
unstrc = 'unstructured'
Models ranked by QIC(x)
15.6.2 Final “Best” Model
<IP_ADDRESS> Tabulate Model Parameters
::knitreg(list(btb_gee_ex),
texregcustom.model.names = c("GEE-ex"),
single.row = TRUE,
caption = "Final GEE, exchangeable correlation")
| GEE-ex | |
|---|---|
| (Intercept) | 5.89 (2.11)** |
| bdi.pre | 0.64 (0.08)*** |
| length>6m | 0.21 (1.48) |
| drugYes | -2.77 (1.65) |
| treatmentBtheB | -2.34 (1.66) |
| month | -0.71 (0.15)*** |
| Scale | 77.14 |
| Num. obs. | 280 |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
<IP_ADDRESS> Plot the Marginal Effects
::interact_plot(model = btb_geeglm_ex,
interactionspred = "month",
modx = "treatment")
Do not worry about confidence intervals.
expand.grid(bdi.pre = 23,
length = "<6m",
drug = "No",
treatment = levels(btb_long$treatment),
month = seq(from = 2, to = 8, by = 2)) %>%
::mutate(fit_in = predict(btb_geeglm_in,
dplyrnewdata = .,
type = "response")) %>%
::mutate(fit_ex = predict(btb_geeglm_ex,
dplyrnewdata = .,
type = "response")) %>%
::mutate(fit_ar = predict(btb_geeglm_ar,
dplyrnewdata = .,
type = "response")) %>%
::mutate(fit_un = predict(btb_geeglm_un,
dplyrnewdata = .,
type = "response")) %>%
::pivot_longer(cols = starts_with("fit_"),
tidyrnames_to = "covR",
names_pattern = "fit_(.*)",
names_ptype = list(covR = "factor()"),
values_to = "fit") %>%
::mutate(covR = factor(covR,
dplyrlevels = c("un", "ar", "ex", "in"),
labels = c("Unstructured",
"Auto-Regressive",
"Compound Symetry",
"Independence"))) %>%
ggplot(aes(x = month,
y = fit,
linetype = treatment)) +
geom_line(alpha = 0.6) +
theme_bw() +
labs(title = "BtheB - Predictions from four GEE models (geeglm)",
x = "Month",
y = "Predicted BDI",
linetype = "Treatment:") +
scale_linetype_manual(values = c("solid", "longdash")) +
theme(legend.key.width = unit(1, "cm")) +
facet_wrap(~ covR)
expand.grid(bdi.pre = 23,
length = "<6m",
drug = "No",
treatment = levels(btb_long$treatment),
month = seq(from = 2, to = 8, by = 2)) %>%
::mutate(fit_in = predict(btb_geeglm_in,
dplyrnewdata = .,
type = "response")) %>%
::mutate(fit_ex = predict(btb_geeglm_ex,
dplyrnewdata = .,
type = "response")) %>%
::mutate(fit_ar = predict(btb_geeglm_ar,
dplyrnewdata = .,
type = "response")) %>%
::mutate(fit_un = predict(btb_geeglm_un,
dplyrnewdata = .,
type = "response")) %>%
::pivot_longer(cols = starts_with("fit_"),
tidyrnames_to = "covR",
names_pattern = "fit_(.*)",
names_ptype = list(covR = "factor()"),
values_to = "fit") %>%
::mutate(covR = factor(covR,
dplyrlevels = c("un", "ar", "ex", "in"),
labels = c("Unstructured",
"Auto-Regressive",
"Compound Symetry",
"Independence"))) %>%
ggplot(aes(x = month,
y = fit,
color = covR,
linetype = treatment)) +
geom_line(alpha = 0.6) +
theme_bw() +
labs(title = "BtheB - Predictions from four GEE models (geeglm)",
x = "Month",
y = "Predicted BDI",
color = "Covariance Pattern:",
linetype = "Treatment:") +
scale_linetype_manual(values = c("solid", "longdash")) +
scale_size_manual(values = c(2, 1, 1, 1)) +
scale_color_manual(values = c("red",
"dodgerblue",
"blue",
"darkgreen")) +
theme(legend.key.width = unit(1, "cm"))
15.6.3 Investigate interactions NOT with time (month
)
<- geepack::geeglm(bdi ~ bdi.pre*length + drug + treatment + month,
btb_geeglm_ex_1 data = btb_long,
id = id,
wave = month,
family = gaussian,
corstr = 'exchangeable')
<- geepack::geeglm(bdi ~ bdi.pre*drug + length + treatment + month,
btb_geeglm_ex_2 data = btb_long,
id = id,
wave = month,
family = gaussian,
corstr = 'exchangeable')
<- geepack::geeglm(bdi ~ bdi.pre*treatment + length + drug + month,
btb_geeglm_ex_3 data = btb_long,
id = id,
wave = month,
family = gaussian,
corstr = 'exchangeable')
::knitreg(list(btb_geeglm_ex,
texreg
btb_geeglm_ex_1,
btb_geeglm_ex_2,
btb_geeglm_ex_3),custom.model.names = c("None",
"Length",
"Antidpressant",
"Treatment"),
label = "GEE_inter1",
caption = "GEE (exchangable): Interactions, not with Time")
| None | Length | Antidpressant | Treatment | |
|---|---|---|---|---|
| (Intercept) | 5.88** | 3.24 | 3.37 | 3.62 |
| (2.11) | (1.99) | (2.68) | (2.78) | |
| bdi.pre | 0.64*** | 0.76*** | 0.75*** | 0.74*** |
| (0.08) | (0.08) | (0.12) | (0.12) | |
| length>6m | 0.20 | 5.66 | 0.47 | 0.22 |
| (1.48) | (3.27) | (1.46) | (1.48) | |
| drugYes | -2.77 | -2.25 | 1.47 | -2.76 |
| (1.65) | (1.62) | (3.20) | (1.64) | |
| treatmentBtheB | -2.33 | -2.42 | -2.19 | 1.24 |
| (1.66) | (1.63) | (1.67) | (3.38) | |
| month | -0.71*** | -0.70*** | -0.71*** | -0.71*** |
| (0.15) | (0.15) | (0.15) | (0.15) | |
| bdi.pre:length>6m | -0.23 | |||
| (0.16) | ||||
| bdi.pre:drugYes | -0.19 | |||
| (0.17) | ||||
| bdi.pre:treatmentBtheB | -0.15 | |||
| (0.17) | ||||
| Scale parameter: gamma | 75.50 | 74.56 | 73.92 | 74.15 |
| Scale parameter: SE | 10.68 | 10.71 | 10.15 | 10.72 |
| Correlation parameter: alpha | 0.69 | 0.69 | 0.68 | 0.68 |
| Correlation parameter: SE | 0.11 | 0.12 | 0.10 | 0.11 |
| Num. obs. | 280 | 280 | 280 | 280 |
| Num. clust. | 97 | 97 | 97 | 97 |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
15.6.4 Investigate interactions with time (month
)
<- geepack::geeglm(bdi ~ bdi.pre + length + drug + treatment*month,
btb_geeglm_ex_11 data = btb_long,
id = id,
wave = month,
family = gaussian,
corstr = 'exchangeable')
<- geepack::geeglm(bdi ~ bdi.pre + length + treatment + drug*month,
btb_geeglm_ex_12 data = btb_long,
id = id,
wave = month,
family = gaussian,
corstr = 'exchangeable')
<- geepack::geeglm(bdi ~ bdi.pre + drug + treatment + length*month,
btb_geeglm_ex_13 data = btb_long,
id = id,
wave = month,
family = gaussian,
corstr = 'exchangeable')
<- geepack::geeglm(bdi ~ length + drug + treatment + bdi.pre*month,
btb_geeglm_ex_14 data = btb_long,
id = id,
wave = month,
family = gaussian,
corstr = 'exchangeable')
::knitreg(list(btb_geeglm_ex,
texreg
btb_geeglm_ex_11,
btb_geeglm_ex_12,
btb_geeglm_ex_13,
btb_geeglm_ex_14),custom.model.names = c("None",
"Treatment",
"Antidepressant",
"Length",
"BL BDI"),
label="GEEinter2",
caption = "GEE (exchangable): Interactions with Time")
| None | Treatment | Antidepressant | Length | BL BDI | |
|---|---|---|---|---|---|
| (Intercept) | 5.88** | 6.83** | 6.63** | 6.28** | 5.09* |
| (2.11) | (2.22) | (2.21) | (2.15) | (2.34) | |
| bdi.pre | 0.64*** | 0.64*** | 0.64*** | 0.64*** | 0.67*** |
| (0.08) | (0.08) | (0.08) | (0.08) | (0.09) | |
| length>6m | 0.20 | 0.25 | 0.16 | -0.44 | 0.27 |
| (1.48) | (1.49) | (1.49) | (1.83) | (1.48) | |
| drugYes | -2.77 | -2.74 | -4.31* | -2.79 | -2.72 |
| (1.65) | (1.66) | (2.01) | (1.64) | (1.66) | |
| treatmentBtheB | -2.33 | -4.23* | -2.32 | -2.31 | -2.36 |
| (1.66) | (2.00) | (1.67) | (1.66) | (1.66) | |
| month | -0.71*** | -0.95*** | -0.88*** | -0.80*** | -0.50 |
| (0.15) | (0.26) | (0.21) | (0.16) | (0.35) | |
| treatmentBtheB:month | 0.47 | ||||
| (0.30) | |||||
| drugYes:month | 0.38 | ||||
| (0.30) | |||||
| length>6m:month | 0.16 | ||||
| (0.29) | |||||
| bdi.pre:month | -0.01 | ||||
| (0.01) | |||||
| Scale parameter: gamma | 75.50 | 76.03 | 76.13 | 75.21 | 75.17 |
| Scale parameter: SE | 10.68 | 10.83 | 10.81 | 10.63 | 10.79 |
| Correlation parameter: alpha | 0.69 | 0.70 | 0.70 | 0.69 | 0.69 |
| Correlation parameter: SE | 0.11 | 0.11 | 0.11 | 0.11 | 0.11 |
| Num. obs. | 280 | 280 | 280 | 280 | 280 |
| Num. clust. | 97 | 97 | 97 | 97 | 97 |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
Now only plot the significant variables for the ‘best’ model.
::interact_plot(model = btb_geeglm_ex,
interactionspred = "month",
modx = "bdi.pre")
::interact_plot(model = btb_geeglm_ex,
interactionspred = "month",
modx = "bdi.pre",
modx.values = c(10, 20, 30),
at = list(length = "<6m",
drug = "No",
treatment = "TAU"),
colors = rep("black", times = 3),
x.label = "Month",
y.label = "Predicted BDI",
legend.main = "Baseline BDI:") +
theme_bw() +
theme(legend.position = "bottom",
legend.key.width = unit(2, "cm")) +
labs(title = "BtheB - Predictions from the GEE model (exchangable)",
subtitle = "Trajectory for a person with BL depression < 6 months and randomized to TAU")
15.7 Conclusion
The Research Question
Does the Beat-the-Blues program
Net of other factors (use of antidepressants and length of the current episode), does the Beat-the-Blues program results in better depression trajectories over treatment as usual?
The Conclusion
There is no evidence that depression trajectories differ between participants randomized to the Beat the Blues program or the Treatment as Usual condition, after accounting for covariates and the correlation between repeated measurements.


This page shows background calculations if the Statistics option in selected in the Background Processing Page of the Options dialog.
CRC32
Calculates a CRC (using standard CRC32 settings) of the whole file. This is the same CRC used on files stored in ZIP format and many other places.
MD5
Calculates the MD5 message digest of the whole file. This format is the most common value currently used to generate message digests and posted on the Internet. However, it is deprecated in favour of SHA1 (below).
SHA1
Calculates the SHA1 message digest of the whole file. Currently the best commonly available message digest. Message digests are most commonly used to verify that files downloaded from public places have not been tampered with.


ABOUT
My name is Cristobal Aguirre. I studied Economics in Chile, worked for 4 years as a senior investment analyst in the country's largest pension fund and then moved to Vancouver with my wife in mid 2017 after 5 months travelling through South East Asia.


public static final java.lang.String |
APPLICATION_DESCRIPTOR_SERVICE_DESCRIPTOR |
"service-descriptor" |
public static final java.lang.String |
CONNECTION_HTTP |
"http" |
public static final java.lang.String |
CONNECTION_HTTPS |
"https" |
public static final java.lang.String |
CONNECTION_QUERY_PARAMETER_KEYWORD |
"?" |
public static final int |
CONNECTION_STATUS_CODE_100 |
100 |
public static final int |
CONNECTION_STATUS_CODE_101 |
101 |
public static final int |
CONNECTION_STATUS_CODE_102 |
102 |
public static final int |
CONNECTION_STATUS_CODE_103 |
103 |
public static final int |
CONNECTION_STATUS_CODE_122 |
122 |
public static final int |
CONNECTION_STATUS_CODE_200 |
200 |
public static final int |
CONNECTION_STATUS_CODE_201 |
201 |
public static final int |
CONNECTION_STATUS_CODE_202 |
202 |
public static final int |
CONNECTION_STATUS_CODE_203 |
203 |
public static final int |
CONNECTION_STATUS_CODE_204 |
204 |
public static final int |
CONNECTION_STATUS_CODE_205 |
205 |
public static final int |
CONNECTION_STATUS_CODE_206 |
206 |
public static final int |
CONNECTION_STATUS_CODE_207 |
207 |
public static final int |
CONNECTION_STATUS_CODE_208 |
208 |
public static final int |
CONNECTION_STATUS_CODE_226 |
226 |
public static final int |
CONNECTION_STATUS_CODE_300 |
300 |
public static final int |
CONNECTION_STATUS_CODE_301 |
301 |
public static final int |
CONNECTION_STATUS_CODE_302 |
302 |
public static final int |
CONNECTION_STATUS_CODE_303 |
303 |
public static final int |
CONNECTION_STATUS_CODE_304 |
304 |
public static final int |
CONNECTION_STATUS_CODE_305 |
305 |
public static final int |
CONNECTION_STATUS_CODE_306 |
306 |
public static final int |
CONNECTION_STATUS_CODE_307 |
307 |
public static final int |
CONNECTION_STATUS_CODE_308 |
308 |
public static final int |
CONNECTION_STATUS_CODE_400 |
400 |
public static final int |
CONNECTION_STATUS_CODE_401 |
401 |
public static final int |
CONNECTION_STATUS_CODE_402 |
402 |
public static final int |
CONNECTION_STATUS_CODE_403 |
403 |
public static final int |
CONNECTION_STATUS_CODE_404 |
404 |
public static final int |
CONNECTION_STATUS_CODE_405 |
405 |
public static final int |
CONNECTION_STATUS_CODE_406 |
406 |
public static final int |
CONNECTION_STATUS_CODE_407 |
407 |
public static final int |
CONNECTION_STATUS_CODE_408 |
408 |
public static final int |
CONNECTION_STATUS_CODE_409 |
409 |
public static final int |
CONNECTION_STATUS_CODE_410 |
410 |
public static final int |
CONNECTION_STATUS_CODE_411 |
411 |
public static final int |
CONNECTION_STATUS_CODE_412 |
412 |
public static final int |
CONNECTION_STATUS_CODE_413 |
413 |
public static final int |
CONNECTION_STATUS_CODE_414 |
414 |
public static final int |
CONNECTION_STATUS_CODE_415 |
415 |
public static final int |
CONNECTION_STATUS_CODE_416 |
416 |
public static final int |
CONNECTION_STATUS_CODE_417 |
417 |
public static final int |
CONNECTION_STATUS_CODE_418 |
418 |
public static final int |
CONNECTION_STATUS_CODE_420 |
420 |
public static final int |
CONNECTION_STATUS_CODE_422 |
422 |
public static final int |
CONNECTION_STATUS_CODE_423 |
423 |
public static final int |
CONNECTION_STATUS_CODE_424 |
424 |
public static final int |
CONNECTION_STATUS_CODE_426 |
426 |
public static final int |
CONNECTION_STATUS_CODE_428 |
428 |
public static final int |
CONNECTION_STATUS_CODE_429 |
429 |
public static final int |
CONNECTION_STATUS_CODE_431 |
431 |
public static final int |
CONNECTION_STATUS_CODE_444 |
444 |
public static final int |
CONNECTION_STATUS_CODE_449 |
449 |
public static final int |
CONNECTION_STATUS_CODE_450 |
450 |
public static final int |
CONNECTION_STATUS_CODE_451 |
451 |
public static final int |
CONNECTION_STATUS_CODE_499 |
499 |
public static final int |
CONNECTION_STATUS_CODE_500 |
500 |
public static final int |
CONNECTION_STATUS_CODE_501 |
501 |
public static final int |
CONNECTION_STATUS_CODE_502 |
502 |
public static final int |
CONNECTION_STATUS_CODE_503 |
503 |
public static final int |
CONNECTION_STATUS_CODE_504 |
504 |
public static final int |
CONNECTION_STATUS_CODE_505 |
505 |
public static final int |
CONNECTION_STATUS_CODE_506 |
506 |
public static final int |
CONNECTION_STATUS_CODE_507 |
507 |
public static final int |
CONNECTION_STATUS_CODE_508 |
508 |
public static final int |
CONNECTION_STATUS_CODE_509 |
509 |
public static final int |
CONNECTION_STATUS_CODE_510 |
510 |
public static final int |
CONNECTION_STATUS_CODE_511 |
511 |
public static final int |
CONNECTION_STATUS_CODE_598 |
598 |
public static final int |
CONNECTION_STATUS_CODE_599 |
599 |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_100 |
"Client should continue with request." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_101 |
"Server is switching protocols." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_102 |
"Server has received and is processing the request." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_103 |
"Resume aborted PUT or POST requests." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_122 |
"URI is longer than a maximum of 2083 characters." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_200 |
"Standard response for successful HTTP requests." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_201 |
"Request has been fulfilled; new resource created." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_202 |
"Request accepted, processing pending." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_203 |
"Request processed, information may be from another source." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_204 |
"Request processed, no content returned." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_205 |
"Request processed, no content returned, reset document view." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_206 |
"Partial resource return due to request header." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_207 |
"XMLl, can contain multiple separate responses." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_208 |
"Results previously returned." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_226 |
"Request fulfilled, response is instance-manipulations." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_300 |
"Multiple options for the resource delivered." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_301 |
"This and all future requests directed to the given URI." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_302 |
"Temporary response to request found via alternative URI." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_303 |
"Permanent response to request found via alternative URI." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_304 |
"Resource has not been modified since last requested." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_305 |
"Content located elsewhere, retrieve from there." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_306 |
"Subsequent requests should use the specified proxy." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_307 |
"Connect again to different URI as provided." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_308 |
"Resumable HTTP Requests." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_400 |
"Request cannot be fulfilled due to bad syntax." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_401 |
"Authentication is possible but has failed." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_402 |
"Payment required, reserved for future use." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_403 |
"Server refuses to respond to request." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_404 |
"Requested resource could not be found." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_405 |
"Request method not supported by that resource." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_406 |
"Content not acceptable according to the Accept headers." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_407 |
"Client must first authenticate itself with the proxy." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_408 |
"Server timed out waiting for the request." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_409 |
"Request could not be processed because of conflict." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_410 |
"Resource is no longer available and will not be available again." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_411 |
"Request did not specify the length of its content." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_412 |
"Server does not meet request preconditions." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_413 |
"Request is larger than the server is willing or able to process." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_414 |
"URI provided was too long for the server to process." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_415 |
"Server does not support media type." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_416 |
"Client has asked for unprovidable portion of the file." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_417 |
"Server cannot meet requirements of Expect request." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_418 |
"I\'m a teapot." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_420 |
"Twitter rate limiting." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_422 |
"Request unable to be followed due to semantic errors." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_423 |
"Resource that is being accessed is locked." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_424 |
"Request failed due to failure of a previous request." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_426 |
"Client should switch to a different protocol." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_428 |
"Origin server requires the request to be conditional." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_429 |
"User has sent too many requests in a given amount of." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_431 |
"Server is unwilling to process the request." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_444 |
"Server returns no information and closes the connection." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_449 |
"Request should be retried after performing action." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_450 |
"Windows Parental Controls blocking access to webpage." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_451 |
"The server cannot reach the client\'s mailbox." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_499 |
"Connection closed by client while HTTP server is processing." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_500 |
"Generic error message." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_501 |
"Server does not recognise method or lacks ability to fulfill." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_502 |
"Server received an invalid response from upstream server." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_503 |
"Server is currently unavailable." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_504 |
"Gateway did not receive response from upstream server." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_505 |
"Server does not support the HTTP protocol version." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_506 |
"Content negotiation for the request results in a circular." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_507 |
"Server is unable to store the representation." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_508 |
"Server detected an infinite loop while processing the request." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_509 |
"Bandwidth limit exceeded." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_510 |
"Further extensions to the request are required." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_511 |
"Client needs to authenticate to gain network access." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_598 |
"Network read timeout behind the proxy." |
public static final java.lang.String |
CONNECTION_STATUS_MESSAGE_599 |
"Network connect timeout behind the proxy." |
public static final java.lang.String |
DATABASE_DESSCRIPTOR_PATH |
"siminov/connect/resource/DatabaseDescriptor.xml" |
public static final java.lang.String |
LIBRARY_DESCRIPTOR_FILE_NAME |
"LibraryDescriptor.xml" |
public static final java.lang.String |
LIBRARY_DESCRIPTOR_PROPERTY |
"property" |
public static final java.lang.String |
LIBRARY_DESCRIPTOR_PROPERTY_NAME |
"name" |
public static final java.lang.String |
LIBRARY_DESCRIPTOR_SERVICE_DESCRIPTOR |
"service-descriptor" |
public static final java.lang.String |
NOTIFICATION_DESCRIPTOR |
"notification-descriptor" |
public static final java.lang.String |
RESOURCE_CLOSE_CURLY_BRACKET |
"}" |
public static final java.lang.String |
RESOURCE_COMMA |
"," |
public static final java.lang.String |
RESOURCE_DOT |
"." |
public static final java.lang.String |
RESOURCE_OPEN_CURLY_BRACKET |
"{" |
public static final java.lang.String |
RESOURCE_REFER_REFERENCE |
"@refer" |
public static final java.lang.String |
RESOURCE_REFERENCE |
"@resource" |
public static final java.lang.String |
RESOURCE_SELF_REFERENCE |
"@self" |
public static final java.lang.String |
RESOURCE_SLASH |
"-" |
public static final java.lang.String |
RESOURCE_SPACE |
" " |
public static final java.lang.String |
SERVICE_DESCRIPTOR_CONTEXT |
"context" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_DESCRIPTION |
"description" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_HTTP_PROTOCOL |
"http" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_HTTPS_PROTOCOL |
"https" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_INSTANCE |
"instance" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_NAME |
"name" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_PORT |
"port" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_PROPERTY |
"property" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_PROPERTY_NAME |
"name" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_PROTOCOL |
"protocol" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST |
"request" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_API |
"api" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_ASYNC_REQUEST_MODE |
"ASYNC" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_CONNECT_TYPE |
"CONNECT" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_DATA_STREAM |
"data-stream" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_DELETE_TYPE |
"DELETE" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_GET_TYPE |
"GET" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_HANDLER |
"handler" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_HEAD_TYPE |
"HEAD" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_HEADER_PARAMETER |
"header-parameter" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_HEADER_PARAMETER_NAME |
"name" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_HEADER_PARAMETER_VALUE |
"value" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_MODE |
"mode" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_NAME |
"name" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_OPTIONS_TYPE |
"OPTIONS" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_PATCH_TYPE |
"PATCH" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_POST_TYPE |
"POST" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_PUT_TYPE |
"PUT" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_QUERY_PARAMETER |
"query-parameter" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_QUERY_PARAMETER_NAME |
"name" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_QUERY_PARAMETER_VALUE |
"value" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_SYNC_REQUEST_MODE |
"SYNC" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_TRACE_TYPE |
"TRACE" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUEST_TYPE |
"type" |
public static final java.lang.String |
SERVICE_DESCRIPTOR_REQUESTS |
"requests" |
public static final java.lang.String |
SYNC_DESCRIPTOR |
"sync-descriptor" |
public static final java.lang.String |
SYNC_DESCRIPTOR_NAME |
"name" |
public static final java.lang.String |
SYNC_DESCRIPTOR_PROPERTY |
"property" |
public static final java.lang.String |
SYNC_DESCRIPTOR_REFRESH_INTERVAL |
"interval" |
public static final java.lang.String |
SYNC_DESCRIPTOR_SERVICE_DESCRIPTOR |
"service-descriptor" |
public static final java.lang.String |
SYNC_DESCRIPTOR_SERVICE_SEPARATOR |
"." |


Denoising Autoencoders (dA)¶
Note
This section assumes the reader has already read through Classifying MNIST digits using Logistic Regression and Multilayer Perceptron. Additionally it uses the following Theano functions and concepts: T.tanh, shared variables, basic arithmetic ops, T.grad, Random numbers, floatX. If you intend to run the code on GPU also read GPU.
Note
The code for this section is available for download here.
The Denoising Autoencoder (dA) is an extension of a classical autoencoder and it was introduced as a building block for deep networks in [Vincent08]. We will start the tutorial with a short discussion on Autoencoders.
Autoencoders¶
See section 4.6 of [Bengio09] for an overview of auto-encoders. An autoencoder takes an input \mathbf{x} \in [0,1]^d and first maps it (with an encoder) to a hidden representation \mathbf{y} \in [0,1]^{d'} through a deterministic mapping, e.g.:
\mathbf{y} = s(\mathbf{W}\mathbf{x} + \mathbf{b})
Where s is a non-linearity such as the sigmoid. The latent representation \mathbf{y}, or code is then mapped back (with a decoder) into a reconstruction \mathbf{z} of the same shape as \mathbf{x}. The mapping happens through a similar transformation, e.g.:
\mathbf{z} = s(\mathbf{W'}\mathbf{y} + \mathbf{b'})
(Here, the prime symbol does not indicate matrix transposition.) \mathbf{z} should be seen as a prediction of \mathbf{x}, given the code \mathbf{y}. Optionally, the weight matrix \mathbf{W'} of the reverse mapping may be constrained to be the transpose of the forward mapping: \mathbf{W'} = \mathbf{W}^T. This is referred to as tied weights. The parameters of this model (namely \mathbf{W}, \mathbf{b}, \mathbf{b'} and, if one doesn’t use tied weights, also \mathbf{W'}) are optimized such that the average reconstruction error is minimized.
The reconstruction error can be measured in many ways, depending on the appropriate distributional assumptions on the input given the code. The traditional squared error L(\mathbf{x} \mathbf{z}) = || \mathbf{x} - \mathbf{z} ||^2, can be used. If the input is interpreted as either bit vectors or vectors of bit probabilities, cross-entropy of the reconstruction can be used:
L_{H} (\mathbf{x}, \mathbf{z}) = - \sum^d_{k=1}[\mathbf{x}_k \log \mathbf{z}_k + (1 - \mathbf{x}_k)\log(1 - \mathbf{z}_k)]
The hope is that the code \mathbf{y} is a distributed representation that captures the coordinates along the main factors of variation in the data. This is similar to the way the projection on principal components would capture the main factors of variation in the data. Indeed, if there is one linear hidden layer (the code) and the mean squared error criterion is used to train the network, then the k hidden units learn to project the input in the span of the first k principal components of the data. If the hidden layer is non-linear, the auto-encoder behaves differently from PCA, with the ability to capture multi-modal aspects of the input distribution. The departure from PCA becomes even more important when we consider stacking multiple encoders (and their corresponding decoders) when building a deep auto-encoder [Hinton06].
Because \mathbf{y} is viewed as a lossy compression of \mathbf{x}, it cannot be a good (small-loss) compression for all \mathbf{x}. Optimization makes it a good compression for training examples, and hopefully for other inputs as well, but not for arbitrary inputs. That is the sense in which an auto-encoder generalizes: it gives low reconstruction error on test examples from the same distribution as the training examples, but generally high reconstruction error on samples randomly chosen from the input space.
We want to implement an auto-encoder using Theano, in the form of a class, that could be afterwards used in constructing a stacked autoencoder. The first step is to create shared variables for the parameters of the autoencoder \mathbf{W}, \mathbf{b} and \mathbf{b'}. (Since we are using tied weights in this tutorial, \mathbf{W}^T will be used for \mathbf{W'}):
def __init__(
self,
numpy_rng,
theano_rng=None,
input=None,
n_visible=784,
n_hidden=500,
W=None,
bhid=None,
bvis=None
):
"""
Initialize the dA class by specifying the number of visible units (the
dimension d of the input ), the number of hidden units ( the dimension
d' of the latent or hidden space ) and the corruption level. The
constructor also receives symbolic variables for the input, weights and
bias. Such a symbolic variables are useful when, for example the input
is the result of some computations, or when weights are shared between
the dA and an MLP layer. When dealing with SdAs this always happens,
the dA on layer 2 gets as input the output of the dA on layer 1,
and the weights of the dA are used in the second stage of training
to construct an MLP.
:type numpy_rng: numpy.random.RandomState
:param numpy_rng: number random generator used to generate weights
:type theano_rng: theano.tensor.shared_randomstreams.RandomStreams
:param theano_rng: Theano random generator; if None is given one is
generated based on a seed drawn from `rng`
:type input: theano.tensor.TensorType
:param input: a symbolic description of the input or None for
standalone dA
:type n_visible: int
:param n_visible: number of visible units
:type n_hidden: int
:param n_hidden: number of hidden units
:type W: theano.tensor.TensorType
:param W: Theano variable pointing to a set of weights that should be
shared belong the dA and another architecture; if dA should
be standalone set this to None
:type bhid: theano.tensor.TensorType
:param bhid: Theano variable pointing to a set of biases values (for
hidden units) that should be shared belong dA and another
architecture; if dA should be standalone set this to None
:type bvis: theano.tensor.TensorType
:param bvis: Theano variable pointing to a set of biases values (for
visible units) that should be shared belong dA and another
architecture; if dA should be standalone set this to None
"""
self.n_visible = n_visible
self.n_hidden = n_hidden
# create a Theano random generator that gives symbolic random values
if not theano_rng:
theano_rng = RandomStreams(numpy_rng.randint(2 ** 30))
# note : W' was written as `W_prime` and b' as `b_prime`
if not W:
# W is initialized with `initial_W` which is uniformely sampled
# from -4*sqrt(6./(n_visible+n_hidden)) and
# 4*sqrt(6./(n_hidden+n_visible))the output of uniform if
# converted using asarray to dtype
# theano.config.floatX so that the code is runable on GPU
initial_W = numpy.asarray(
numpy_rng.uniform(
low=-4 * numpy.sqrt(6. / (n_hidden + n_visible)),
high=4 * numpy.sqrt(6. / (n_hidden + n_visible)),
size=(n_visible, n_hidden)
),
dtype=theano.config.floatX
)
W = theano.shared(value=initial_W, name='W', borrow=True)
if not bvis:
bvis = theano.shared(
value=numpy.zeros(
n_visible,
dtype=theano.config.floatX
),
borrow=True
)
if not bhid:
bhid = theano.shared(
value=numpy.zeros(
n_hidden,
dtype=theano.config.floatX
),
name='b',
borrow=True
)
self.W = W
# b corresponds to the bias of the hidden
self.b = bhid
# b_prime corresponds to the bias of the visible
self.b_prime = bvis
# tied weights, therefore W_prime is W transpose
self.W_prime = self.W.T
self.theano_rng = theano_rng
# if no input is given, generate a variable representing the input
if input is None:
# we use a matrix because we expect a minibatch of several
# examples, each example being a row
self.x = T.dmatrix(name='input')
else:
self.x = input
self.params = [self.W, self.b, self.b_prime]
Note that we pass the symbolic input
to the autoencoder as a parameter.
This is so that we can concatenate layers of autoencoders to form a deep
network: the symbolic output (the \mathbf{y} above) of layer k will
be the symbolic input of layer k+1.
Now we can express the computation of the latent representation and of the reconstructed signal:
def get_hidden_values(self, input):
""" Computes the values of the hidden layer """
return T.nnet.sigmoid(T.dot(input, self.W) + self.b)
def get_reconstructed_input(self, hidden):
"""Computes the reconstructed input given the values of the
hidden layer
"""
return T.nnet.sigmoid(T.dot(hidden, self.W_prime) + self.b_prime)
And using these functions we can compute the cost and the updates of one stochastic gradient descent step:
def get_cost_updates(self, corruption_level, learning_rate):
""" This function computes the cost and the updates for one trainng
step of the dA """
tilde_x = self.get_corrupted_input(self.x, corruption_level)
y = self.get_hidden_values(tilde_x)
z = self.get_reconstructed_input(y)
# note : we sum over the size of a datapoint; if we are using
# minibatches, L will be a vector, with one entry per
# example in minibatch
L = - T.sum(self.x * T.log(z) + (1 - self.x) * T.log(1 - z), axis=1)
# note : L is now a vector, where each element is the
# cross-entropy cost of the reconstruction of the
# corresponding example of the minibatch. We need to
# compute the average of all these to get the cost of
# the minibatch
cost = T.mean(L)
# compute the gradients of the cost of the `dA` with respect
# to its parameters
gparams = T.grad(cost, self.params)
# generate the list of updates
updates = [
(param, param - learning_rate * gparam)
for param, gparam in zip(self.params, gparams)
]
return (cost, updates)
We can now define a function that applied iteratively will update the
parameters W
, b
and b_prime
such that the
reconstruction cost is approximately minimized.
da = dA(
numpy_rng=rng,
theano_rng=theano_rng,
input=x,
n_visible=28 * 28,
n_hidden=500
)
cost, updates = da.get_cost_updates(
corruption_level=0.,
learning_rate=learning_rate
)
train_da = theano.function(
[index],
cost,
updates=updates,
givens={
x: train_set_x[index * batch_size: (index + 1) * batch_size]
}
)
start_time = timeit.default_timer()
############
# TRAINING #
############
# go through training epochs
for epoch in range(training_epochs):
# go through trainng set
c = []
for batch_index in range(n_train_batches):
c.append(train_da(batch_index))
print('Training epoch %d, cost ' % epoch, numpy.mean(c, dtype='float64'))
end_time = timeit.default_timer()
training_time = (end_time - start_time)
print(('The no corruption code for file ' +
os.path.split(__file__)[1] +
' ran for %.2fm' % ((training_time) / 60.)), file=sys.stderr)
image = Image.fromarray(
tile_raster_images(X=da.W.get_value(borrow=True).T,
img_shape=(28, 28), tile_shape=(10, 10),
tile_spacing=(1, 1)))
image.save('filters_corruption_0.png')
# start-snippet-3
#####################################
# BUILDING THE MODEL CORRUPTION 30% #
#####################################
rng = numpy.random.RandomState(123)
theano_rng = RandomStreams(rng.randint(2 ** 30))
da = dA(
numpy_rng=rng,
theano_rng=theano_rng,
input=x,
n_visible=28 * 28,
n_hidden=500
)
cost, updates = da.get_cost_updates(
corruption_level=0.3,
learning_rate=learning_rate
)
train_da = theano.function(
[index],
cost,
updates=updates,
givens={
x: train_set_x[index * batch_size: (index + 1) * batch_size]
}
)
start_time = timeit.default_timer()
############
# TRAINING #
############
# go through training epochs
for epoch in range(training_epochs):
# go through trainng set
c = []
for batch_index in range(n_train_batches):
c.append(train_da(batch_index))
print('Training epoch %d, cost ' % epoch, numpy.mean(c, dtype='float64'))
end_time = timeit.default_timer()
training_time = (end_time - start_time)
print(('The 30% corruption code for file ' +
os.path.split(__file__)[1] +
' ran for %.2fm' % (training_time / 60.)), file=sys.stderr)
# end-snippet-3
# start-snippet-4
image = Image.fromarray(tile_raster_images(
X=da.W.get_value(borrow=True).T,
img_shape=(28, 28), tile_shape=(10, 10),
tile_spacing=(1, 1)))
image.save('filters_corruption_30.png')
# end-snippet-4
os.chdir('../')
if __name__ == '__main__':
test_dA()
If there is no constraint besides minimizing the reconstruction error, one might expect an auto-encoder with n inputs and an encoding of dimension n (or greater) to learn the identity function, merely mapping an input to its copy. Such an autoencoder would not differentiate test examples (from the training distribution) from other input configurations.
Surprisingly, experiments reported in [Bengio07] suggest that, in practice, when trained with stochastic gradient descent, non-linear auto-encoders with more hidden units than inputs (called overcomplete) yield useful representations. (Here, “useful” means that a network taking the encoding as input has low classification error.)
A simple explanation is that stochastic gradient descent with early stopping is similar to an L2 regularization of the parameters. To achieve perfect reconstruction of continuous inputs, a one-hidden layer auto-encoder with non-linear hidden units (exactly like in the above code) needs very small weights in the first (encoding) layer, to bring the non-linearity of the hidden units into their linear regime, and very large weights in the second (decoding) layer. With binary inputs, very large weights are also needed to completely minimize the reconstruction error. Since the implicit or explicit regularization makes it difficult to reach large-weight solutions, the optimization algorithm finds encodings which only work well for examples similar to those in the training set, which is what we want. It means that the representation is exploiting statistical regularities present in the training set, rather than merely learning to replicate the input.
There are other ways by which an auto-encoder with more hidden units than inputs could be prevented from learning the identity function, capturing something useful about the input in its hidden representation. One is the addition of sparsity (forcing many of the hidden units to be zero or near-zero). Sparsity has been exploited very successfully by many [Ranzato07] [Lee08]. Another is to add randomness in the transformation from input to reconstruction. This technique is used in Restricted Boltzmann Machines (discussed later in Restricted Boltzmann Machines (RBM)), as well as in Denoising Auto-Encoders, discussed below.
Denoising Autoencoders¶
The idea behind denoising autoencoders is simple. In order to force the hidden layer to discover more robust features and prevent it from simply learning the identity, we train the autoencoder to reconstruct the input from a corrupted version of it.
The denoising auto-encoder is a stochastic version of the auto-encoder. Intuitively, a denoising auto-encoder does two things: try to encode the input (preserve the information about the input), and try to undo the effect of a corruption process stochastically applied to the input of the auto-encoder. The latter can only be done by capturing the statistical dependencies between the inputs. The denoising auto-encoder can be understood from different perspectives (the manifold learning perspective, stochastic operator perspective, bottom-up – information theoretic perspective, top-down – generative model perspective), all of which are explained in [Vincent08]. See also section 7.2 of [Bengio09] for an overview of auto-encoders.
In [Vincent08], the stochastic corruption process randomly sets some of the inputs (as many as half of them) to zero. Hence the denoising auto-encoder is trying to predict the corrupted (i.e. missing) values from the uncorrupted (i.e., non-missing) values, for randomly selected subsets of missing patterns. Note how being able to predict any subset of variables from the rest is a sufficient condition for completely capturing the joint distribution between a set of variables (this is how Gibbs sampling works).
To convert the autoencoder class into a denoising autoencoder class, all we need to do is to add a stochastic corruption step operating on the input. The input can be corrupted in many ways, but in this tutorial we will stick to the original corruption mechanism of randomly masking entries of the input by making them zero. The code below does just that:
def get_corrupted_input(self, input, corruption_level):
"""This function keeps ``1-corruption_level`` entries of the inputs the
same and zero-out randomly selected subset of size ``corruption_level``
Note : first argument of theano.rng.binomial is the shape(size) of
random numbers that it should produce
second argument is the number of trials
third argument is the probability of success of any trial
this will produce an array of 0s and 1s where 1 has a
probability of 1 - ``corruption_level`` and 0 with
``corruption_level``
The binomial function return int64 data type by
default. int64 multiplicated by the input
type(floatX) always return float64. To keep all data
in floatX when floatX is float32, we set the dtype of
the binomial to floatX. As in our case the value of
the binomial is always 0 or 1, this don't change the
result. This is needed to allow the gpu to work
correctly as it only support float32 for now.
"""
return self.theano_rng.binomial(size=input.shape, n=1,
p=1 - corruption_level,
dtype=theano.config.floatX) * input
In the stacked autoencoder class (Stacked Autoencoders) the weights of
the dA
class have to be shared with those of a corresponding sigmoid layer.
For this reason, the constructor of the dA
also gets Theano variables
pointing to the shared parameters. If those parameters are left to None
,
new ones will be constructed.
The final denoising autoencoder class becomes:
class dA(object):
"""Denoising Auto-Encoder class (dA)
A denoising autoencoders tries to reconstruct the input from a corrupted
version of it by projecting it first in a latent space and reprojecting
it afterwards back in the input space. Please refer to Vincent et al.,2008
for more details. If x is the input then equation (1) computes a partially
destroyed version of x by means of a stochastic mapping q_D. Equation (2)
computes the projection of the input into the latent space. Equation (3)
computes the reconstruction of the input, while equation (4) computes the
reconstruction error.
.. math::
\tilde{x} ~ q_D(\tilde{x}|x) (1)
y = s(W \tilde{x} + b) (2)
x = s(W' y + b') (3)
L(x,z) = -sum_{k=1}^d [x_k \log z_k + (1-x_k) \log( 1-z_k)] (4)
"""
def __init__(
self,
numpy_rng,
theano_rng=None,
input=None,
n_visible=784,
n_hidden=500,
W=None,
bhid=None,
bvis=None
):
"""
Initialize the dA class by specifying the number of visible units (the
dimension d of the input ), the number of hidden units ( the dimension
d' of the latent or hidden space ) and the corruption level. The
constructor also receives symbolic variables for the input, weights and
bias. Such a symbolic variables are useful when, for example the input
is the result of some computations, or when weights are shared between
the dA and an MLP layer. When dealing with SdAs this always happens,
the dA on layer 2 gets as input the output of the dA on layer 1,
and the weights of the dA are used in the second stage of training
to construct an MLP.
:type numpy_rng: numpy.random.RandomState
:param numpy_rng: number random generator used to generate weights
:type theano_rng: theano.tensor.shared_randomstreams.RandomStreams
:param theano_rng: Theano random generator; if None is given one is
generated based on a seed drawn from `rng`
:type input: theano.tensor.TensorType
:param input: a symbolic description of the input or None for
standalone dA
:type n_visible: int
:param n_visible: number of visible units
:type n_hidden: int
:param n_hidden: number of hidden units
:type W: theano.tensor.TensorType
:param W: Theano variable pointing to a set of weights that should be
shared belong the dA and another architecture; if dA should
be standalone set this to None
:type bhid: theano.tensor.TensorType
:param bhid: Theano variable pointing to a set of biases values (for
hidden units) that should be shared belong dA and another
architecture; if dA should be standalone set this to None
:type bvis: theano.tensor.TensorType
:param bvis: Theano variable pointing to a set of biases values (for
visible units) that should be shared belong dA and another
architecture; if dA should be standalone set this to None
"""
self.n_visible = n_visible
self.n_hidden = n_hidden
# create a Theano random generator that gives symbolic random values
if not theano_rng:
theano_rng = RandomStreams(numpy_rng.randint(2 ** 30))
# note : W' was written as `W_prime` and b' as `b_prime`
if not W:
# W is initialized with `initial_W` which is uniformely sampled
# from -4*sqrt(6./(n_visible+n_hidden)) and
# 4*sqrt(6./(n_hidden+n_visible))the output of uniform if
# converted using asarray to dtype
# theano.config.floatX so that the code is runable on GPU
initial_W = numpy.asarray(
numpy_rng.uniform(
low=-4 * numpy.sqrt(6. / (n_hidden + n_visible)),
high=4 * numpy.sqrt(6. / (n_hidden + n_visible)),
size=(n_visible, n_hidden)
),
dtype=theano.config.floatX
)
W = theano.shared(value=initial_W, name='W', borrow=True)
if not bvis:
bvis = theano.shared(
value=numpy.zeros(
n_visible,
dtype=theano.config.floatX
),
borrow=True
)
if not bhid:
bhid = theano.shared(
value=numpy.zeros(
n_hidden,
dtype=theano.config.floatX
),
name='b',
borrow=True
)
self.W = W
# b corresponds to the bias of the hidden
self.b = bhid
# b_prime corresponds to the bias of the visible
self.b_prime = bvis
# tied weights, therefore W_prime is W transpose
self.W_prime = self.W.T
self.theano_rng = theano_rng
# if no input is given, generate a variable representing the input
if input is None:
# we use a matrix because we expect a minibatch of several
# examples, each example being a row
self.x = T.dmatrix(name='input')
else:
self.x = input
self.params = [self.W, self.b, self.b_prime]
def get_corrupted_input(self, input, corruption_level):
"""This function keeps ``1-corruption_level`` entries of the inputs the
same and zero-out randomly selected subset of size ``corruption_level``
Note : first argument of theano.rng.binomial is the shape(size) of
random numbers that it should produce
second argument is the number of trials
third argument is the probability of success of any trial
this will produce an array of 0s and 1s where 1 has a
probability of 1 - ``corruption_level`` and 0 with
``corruption_level``
The binomial function return int64 data type by
default. int64 multiplicated by the input
type(floatX) always return float64. To keep all data
in floatX when floatX is float32, we set the dtype of
the binomial to floatX. As in our case the value of
the binomial is always 0 or 1, this don't change the
result. This is needed to allow the gpu to work
correctly as it only support float32 for now.
"""
return self.theano_rng.binomial(size=input.shape, n=1,
p=1 - corruption_level,
dtype=theano.config.floatX) * input
def get_hidden_values(self, input):
""" Computes the values of the hidden layer """
return T.nnet.sigmoid(T.dot(input, self.W) + self.b)
def get_reconstructed_input(self, hidden):
"""Computes the reconstructed input given the values of the
hidden layer
"""
return T.nnet.sigmoid(T.dot(hidden, self.W_prime) + self.b_prime)
def get_cost_updates(self, corruption_level, learning_rate):
""" This function computes the cost and the updates for one trainng
step of the dA """
tilde_x = self.get_corrupted_input(self.x, corruption_level)
y = self.get_hidden_values(tilde_x)
z = self.get_reconstructed_input(y)
# note : we sum over the size of a datapoint; if we are using
# minibatches, L will be a vector, with one entry per
# example in minibatch
L = - T.sum(self.x * T.log(z) + (1 - self.x) * T.log(1 - z), axis=1)
# note : L is now a vector, where each element is the
# cross-entropy cost of the reconstruction of the
# corresponding example of the minibatch. We need to
# compute the average of all these to get the cost of
# the minibatch
cost = T.mean(L)
# compute the gradients of the cost of the `dA` with respect
# to its parameters
gparams = T.grad(cost, self.params)
# generate the list of updates
updates = [
(param, param - learning_rate * gparam)
for param, gparam in zip(self.params, gparams)
]
return (cost, updates)
Putting it All Together¶
It is easy now to construct an instance of our dA
class and train
it.
# allocate symbolic variables for the data
index = T.lscalar() # index to a [mini]batch
x = T.matrix('x') # the data is presented as rasterized images
#####################################
# BUILDING THE MODEL CORRUPTION 30% #
#####################################
rng = numpy.random.RandomState(123)
theano_rng = RandomStreams(rng.randint(2 ** 30))
da = dA(
numpy_rng=rng,
theano_rng=theano_rng,
input=x,
n_visible=28 * 28,
n_hidden=500
)
cost, updates = da.get_cost_updates(
corruption_level=0.3,
learning_rate=learning_rate
)
train_da = theano.function(
[index],
cost,
updates=updates,
givens={
x: train_set_x[index * batch_size: (index + 1) * batch_size]
}
)
start_time = timeit.default_timer()
############
# TRAINING #
############
# go through training epochs
for epoch in range(training_epochs):
# go through trainng set
c = []
for batch_index in range(n_train_batches):
c.append(train_da(batch_index))
print('Training epoch %d, cost ' % epoch, numpy.mean(c, dtype='float64'))
end_time = timeit.default_timer()
training_time = (end_time - start_time)
print(('The 30% corruption code for file ' +
os.path.split(__file__)[1] +
' ran for %.2fm' % (training_time / 60.)), file=sys.stderr)
In order to get a feeling of what the network learned we are going to plot the filters (defined by the weight matrix). Bear in mind, however, that this does not provide the entire story, since we neglect the biases and plot the weights up to a multiplicative constant (weights are converted to values between 0 and 1).
To plot our filters we will need the help of tile_raster_images
(see
Plotting Samples and Filters) so we urge the reader to study it. Also
using the help of the Python Image Library, the following lines of code will
save the filters as an image:
image = Image.fromarray(tile_raster_images(
X=da.W.get_value(borrow=True).T,
img_shape=(28, 28), tile_shape=(10, 10),
tile_spacing=(1, 1)))
image.save('filters_corruption_30.png')
Running the Code¶
To run the code:
python dA.py
The resulted filters when we do not use any noise are:
The filters for 30 percent noise:


Suspect arrested in Khanewal for murdering man, burning his body 'in the name of honour'
by https://www.dawn.com/authors/5579/sajjad-akbar-shah, Sajjad Akbar ShahPolice in Jahanian tehsil of Khanewal district on Thursday found the body of a man who was killed and burnt by an individual allegedly in the name of honour, said Deputy Superintendent of Police (DSP) Jahanian Shamsuddin.
According to the DSP, the suspect had tortured and murdered the victim, and then burnt his body.
As per the initial investigation, the DSP said the suspect had suspicions of illicit relations between his cousin and the victim. He added that the murder was allegedly carried out in the name of honour.
Shamsuddin said that the victim's body has been shifted to a hospital while further investigation into the incident is underway.
Meanwhile, Khanewal District Police Officer (DPO) Faisal Shehzad recognised the timely response of DSP Jahanian and his team in arresting the suspect within an hour of receiving information.
According to DPO Shehzad, the police had received information that the burnt body of a man had been found in Jahanian after which the DSP, his Crime Investigation Agency (CIA) and forensic teams responded and cordoned off the site of the incident.
Following this, the forensic team started to collect evidence and began an interrogation with local people in the area during which they received a clue that led to the arrest of the main suspect.
Khanewal DPO said that their investigation was underway in order to determine the motive behind the murder, the weapon used as well as in order to determine who else was involved in the incident.
He added that the suspect had admitted to the crime during the investigation.


Please Contact Us With Any Questions
Fifth Capital Management
3131 McKinney Avenue
Suite 600
Dallas, TX 75204
800.766.4997
Fifth Capital Management
3131 McKinney Avenue
Suite 600
Dallas, TX 75204
800.766.4997


{% extends "layouts/layout1.html" %} {% set active_page = "login_user" %} {% block title %}Login{% endblock %} {% block content %}
Coalman is a collaboration platform that helps manage the bid process for your consturction projects.
Forgot Password


Contains all actions for the engine that can alter different sections of a game including Entities and Levels.


{% extends "base.html" %} {% load staticfiles %} {% block content %}
Welcome, {{ user.name }}This is a mini-blog application using a back-end built on Django 2.0 and Django Rest Framework. It illustrates how to create and send JSON Web Token authentication headers with the HTTP requests.