Plotting the spline model .pull-left[ We cannot use `geom_smooth` to show the fit of this spline model. However, with some additional effort, we can plot the effect on a graph. 1. Calculate the predicted value of wages for a range of age values based on the model using the `predict` command. 2. Feed in these predicted values as a dataset to the `geom_line` function in our `ggplot` graph. ] -- .pull-right[ ### The `predict` command The `predict` command takes two important arguments: 1. The model object for which we want predicted values of the dependent variable. 2. A new dataset that contains all the same independent variables that are in the model. ```r model <- lm(wages~age+age.spline, data=earnings) pre_df <- data.frame(age=18:65) pre_df$age.spline <- ifelse(pre_df$age<35, 0, pre_df$age-35) pre_df$wages <- predict(model, newdata=pre_df) head(pre_df, n=3) ``` ``` ## age age.spline wages ## 1 18 0 11.01100 ## 2 19 0 11.95824 ## 3 20 0 12.90548 ``` ] --- ## Adding the spline fit to ggplot .pull-left[ ```r ggplot(earnings, aes(x=age, y=wages))+ geom_jitter(alpha=0.01, width = 1)+ * geom_line(data=predict_df, color="blue", size=1.5) ``` To add the spline fit, we just add a `geom_line` command and specify our `pre_df` dataset through the `data` argument so that it uses the predicted values rather than the actual `earnings` dataset to graph the line. ] .pull-right[  ] --- class: inverse, center, middle background-image: url(images/bence-balla-schottner-o5ttfkZU9_Q-unsplash.jpg) background-size: cover name: iid # The IID Violation and Robust Standard Errors --- ## Linear model as data-generating process `$$y_i=\underbrace{\beta_0+\beta_1x_{i1}+\beta_2x_{i2}+\ldots+\beta_px_{ip}}_{\hat{y}_i}+\epsilon_i$$` -- - To get values of `\(y_i\)`, you feed in values of `\(x_i\)` to the structural component and get back out a predicted `\(\hat{y_i}\)` value. -- - To get the stochastic (random) part, you then reach into some distribution to grab a random value of `\(\epsilon_i\)` that you add to your predicted value to get an actual value of `\(y_i\)`. -- - What distribution are you reaching into when you grab `\(\epsilon_i\)`? The linear model framework assumes two things about this distribution: -- - **Independence**: The number you pull out each time doesn't depend on other numbers that you pull out. This is most commonly violated by **autocorrelation** or the clustering of **repeated observations**. -- - **Identical distribution**: You reach into the same distribution for all observations. This is most commonly violated by **heteroscedasticity** which means non-constant variance in the residuals. -- - Together these assumptions give us the **IID** assumption of the linear model: the error terms are **independent and identically distributed**. --- ## Violations of independence -- .pull-left[ ### Serial autocorrelation  In time series data, sequential observations in time are likely to either be highly positively or negatively correlated. In this case, when one year has a higher or lower than expected value, the next and previous years are also likely to have a higher or lower than expected value. ] -- .pull-right[ ### Repeated observations .center[![](images/classroom.jpg)] When observations are drawn repeatedly from a sample of some larger units (i.e. a multilevel structure), then observations within the same unit are likely to vary from predicted values in the same way. For example, students in the same classroom might tend to have either lower or higher test scores than predicted by the model due to some unobserved feature of that classroom. ] --- ## Serial autocorrelation example .pull-left[ ```r model <- lm(Employed~GNP, data=longley) temp <- augment(model)$.resid temp <- data.frame(x=temp[1:15], y=temp[2:16]) ggplot(temp, aes(x=x, y=y))+ geom_point()+ geom_smooth(method="lm", se=FALSE)+ geom_hline(yintercept = 0, linetype=2)+ labs(x="residuals from 1947-1961", y="residuals from 1948-1962")+ annotate("text",0.25, 0.5, label="r=0.16") ``` As an example of serial autocorrelation, I will use the `longley` time series dataset in R to fit the following model predicting the number of people employed by GNP from 1947 to 1962 (n=16): I can then plot the residuals values from years 1947 to 1961 by the residual values for years 1948 to 1962. The positive correlation in the residuals here suggests serial autocorrelation. ] .pull-right[  ] --- ## Heteroscedasticity .pull-left[ .center[![](images/mary_poppins.jpg)] - Heteroscedasticity means that the variance of the residuals is not constant but depends on the values of `\(x_i\)`, and therefore, implicitly, `\(\hat{y}_i\)`. - A classic example of heteroscedasticity is when the variance of the residuals increases with larger values of `\(\hat{y}_i\)` giving you a cone shape in a residual by fitted value plot. ] .pull-right[  ] --- ## Correcting for iid violations Violating the iid. assumption does not bias your results, but it will lead to inefficient estimates and poorly estimated standard errors. There are a number of potential solutions to the iid problem. These include: -- - **Transformations** (particularly the log transformation) can often solve the problem of heteroscedasticity. -- - **Weighted least squares** models can correct for iid when the nature of the violation is understood. -- - **Robust standard errors** can be used as a crude brute-force solution when the nature of the violation is not well understood. I would recommend that this only be done for diagnostic reasons. -- - In general, the best approach is to re-think your model. If you have an iid violation then you are probably not applying the best type of model to the problem at hand. --- ## Fixing heteroscedasticity with a transformation .pull-left[  ] -- .pull-right[  ] --- ## Using Weighted Least Squares The weighted least squares technique uses a weighting matrix `\(\mathbf{W}\)` in its calculation of regression slopes like so: `$$\mathbf{\beta}=(\mathbf{X'W^{-1}X})^{-1}\mathbf{X'W^{-1}y}$$` `$$SE_{\beta}=\sqrt{\sigma^{2}(\mathbf{X'W^{-1}X})^{-1}}$$` The exact form of this weighting matrix depends on the nature of the iid violation, but in general it is used to represent the covariance between residuals. Values in the diagonal cells adjust for heteroscedasticity and values in other cells adjust for autocorrelation. This technique is used routinely in time series analysis to adjust for serial autocorrelation. --- ## GLS Example We can use the `gls` command in the `nmle` package to adjust for serial autocorrelation in the prior model.We will assume that the autocorrelation follows an "AR1" pattern in which each subsequent residual is only correlated with its immediate predecessor (AR1 stands for auto-regressive 1, where 1 indicates the lag). ```r summary(lm(Employed~GNP+Population, data=longley))$coef ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 88.93879831 13.78502699 6.451841 2.160219e-05 ## GNP 0.06317244 0.01064721 5.933238 4.958268e-05 ## Population -0.40974292 0.15213679 -2.693253 1.842949e-02 ``` ```r library(nlme) model.ar1 <- gls(Employed~GNP+Population, correlation=corAR1(form=~Year),data=longley) summary(model.ar1)$tTable ``` ``` ## Value Std.Error t-value p-value ## (Intercept)<PHONE_NUMBER>5 14.1989324 7.173647 7.222049e-06 ## GNP 0.07207088 0.0106057 6.795485 1.271171e-05 ## Population -0.54851350 0.1541297 -3.558778 3.497130e-03 ``` --- ## Robust Standard Errors .pull-left[ We won't delve into the math behind the robust standard error, but the general idea is that robust standard errors will give you "correct" standard errors even when the model is mis-specified due to issues such a non-linearity, heteroscedasticity, and autocorrelation. Robust standard errors can be estimated in R using the `sandwich` and `lmtest` packages, and specifically with the `coeftest` command. Within this command, it is possible to specify different types of robust standard errors, but we will use the "HC1" version which is equivalent to the robust standard errors produced in Stata by default. ] -- .pull-right[ ```r library(sandwich) library(lmtest) model <- lm(BoxOffice~TomatoRating, data=movies) summary(model)$coef[,1:2] ``` ``` ## Estimate Std. Error ## (Intercept) -16.71941 5.1103643 ## TomatoRating 11.43015 0.9140331 ``` ```r coeftest(model, vcov=vcovHC(model, "HC1"))[,1:2] ``` ``` ## Estimate Std. Error ## (Intercept) -16.71941 5.257507 ## TomatoRating 11.43015 1.097638 ``` Note that the estimates are the same, but the robust standard errors are considerably larger. That difference in magnitude is telling us that our basic regression model is problematic. In this case, we already know that the problem is heteroscedasticity. ] --- ## ⚠ Robust standard errors diagnose problems that they do not fix! ```r model <- lm(log(BoxOffice)~TomatoRating, data=movies) summary(model)$coef ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.9727622 0.20091999 4.841540 1.365329e-06 ## TomatoRating 0.2406026 0.03593629 6.695254 2.638277e-11 ``` ```r coeftest(model, vcov=vcovHC(model, "HC1")) ``` ``` ## ## t test of coefficients: ## ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.972762 0.186379 5.2193 1.941e-07 *** ## TomatoRating 0.240603 0.032115 7.4919 9.302e-14 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` When the dependent variable is logged to remove heteroscedasticity, the difference between robust and regular standard errors goes away. --- class: inverse, center, middle background-image: url(images/cyril-saulnier-TsVN31Dzyv4-unsplash.jpg) background-size: cover name: sample-design # Sample Design and Weighting --- ## The reality of survey sampling -- .pull-left[ ### The simple random sample .center[![](images/srs.png)] - We have treated all of our sample data as if it was collected using simple random sampling (SRS). - In an SRS of size `\(n\)`, every possible combination of `\(n\)` observations from the population has an equally likely chance of being drawn. - In an SRS, every statistic (e.g. mean, proportions, linear slopes) should be representative of the population, except for random sampling bias. ] -- .pull-right[ ### Reality .center[![](images/strat_sample.png)] - In practice, large-scale surveys never use SRS for pragmatic and design reasons. - For correct statistical inference, we typically need to make adjustments for sample design. - The primary ways in which sample design affects estimation are**clustering**, **stratification**, and **weighting**. ] --- ## Cluster/Multistage sampling .left-column[ ![](images/multistage.png) ] .right-column[ Potential observations are aggregated into larger groupings identified as the Primary Sampling Unit (PSU) and then we sample some of these PSUs before sampling individual observations within the sampled PSUs. - Cluster sampling is about efficiency and cost. Clusters are typically defined geographically, which then minimizes the cost of sampling individual observations. - If PSUs are sampled with probabilities proportional to cluster size, then every unit in the population has an equal likelihood of being selected. In this case, summary statistics on the sample should be representative. - When the variable of interest is distributed differently across clusters, the sampling variability will be higher than an SRS even if every observation has an equally likely chance of being drawn. ] --- ## What percent of the US population is Mormon? .pull-left[ The General Social Survey uses Metropolitan Statistical Areas (MSA) and clusters of non-metropolitan counties as the PSU. In each year, it draws a sample of these areas and then samples respondents within each sampled PSU. Because Mormons are heavily concentrated in certain places, percent Mormon in the GSS varies substantially from year to year by whether certain PSUs were sampled or not. Red band indicates expected 95% interval for sampling variability without clustering, assuming the average percent across all years (dotted line). ] .pull-right[  ] --- ## Stratification Stratification in sampling operates in a manner somewhat similar to cluster sampling except that once the observations are aggregated into strata by some characteristic (e.g. income, race, age), observations are sampled from *every* stratum. -- - Stratification is typically done to ensure that various sub-populations are present within the sample. -- - Different strata may be sampled with different probabilities. The most common approach is to take an **oversample** of a small group in order to ensure that effective comparisons can be made between that group and other groups. -- - In practice, stratification is often done by first screening potential respondents for stratum characteristics. -- - If strata are sampled with different probabilities, then summary statistics for the full sample will not be representative without weight adjustments. -- - Unlike clustering, greater similarity on a characteristic of interest within strata can actually **reduce** the sampling variability for that characteristic relative to an SRS. --- ## Weighting -- - Numerous factors can lead to a sample being unrepresentative such as sampling strata with different probabilities, differential non-response rates, and a lack of fit between sample frame and population. -- - **Sampling weights** allow researchers to correct summary statistics from the sample so that they are representative of the population. -- - The sampling weight for an observation should be `\(1/p_i\)` where `\(p_i\)` is the probability of being sampled. The sampling weight indicates the number of observations in the population that observation `\(i\)` in the sample represents. -- - Calculating sampling weights can be quite complex. In some cases, researchers may know `\(p_i\)` from the study design. In other cases, researchers may create **post-stratification** weights by comparing the sample to some other data source (e.g. census, school records) for a set of demographic characteristics and applying weights to make the sample align with the other data source. -- - When sampling weights are present in a dataset, they must be used to generate statistics representative of the population. -- - Variation in sampling weights will increase sampling variability above and beyond that expected for an SRS. --- ## Weights in Sneetchville .left-column[ ![sneetches](images/sneetches.jpg) ] .right-column[ In Sneetchville, there are 3 star-bellied and 7 regular sneetches. Lets say I take a stratified sample of two star-bellied and two regular sneetches. Here is the population data: ```r sneetches <- data.frame(type=factor(c(rep("Star-bellied", 3),rep("Regular",7))), income=c(3,7,6,2,1,4,5,0,2,1)) sneetches$prob <- c(rep(2/3,3),rep(2/7,7)) sneetches ``` ``` ## type income prob ## 1 Star-bellied 3 0.6666667 ## 2 Star-bellied 7 0.6666667 ## 3 Star-bellied 6 0.6666667 ## 4 Regular 2 0.2857143 ## 5 Regular 1 0.2857143 ## 6 Regular 4 0.2857143 ## 7 Regular 5 0.2857143 ## 8 Regular 0 0.2857143 ## 9 Regular 2 0.2857143 ## 10 Regular 1 0.2857143 ``` ] --- ## Weights in Sneetchville ```r sneetch.sample <- sneetches[c(sample(1:3,2),sample(4:10,2)),] sneetch.sample$weight <- 1/sneetch.sample$prob sneetch.sample ``` ``` ## type income prob weight ## 2 Star-bellied 7 0.6666667 1.5 ## 3 Star-bellied 6 0.6666667 1.5 ## 7 Regular 5 0.2857143 3.5 ## 6 Regular 4 0.2857143 3.5 ``` ```r x <- c(mean(sneetches$income), mean(sneetch.sample$income), weighted.mean(sneetch.sample$income, sneetch.sample$weight)) names(x) <- c("population","sample.unweighted","sample.weighted") x ``` ``` ## population sample.unweighted sample.weighted ## 3.1 5.5 5.1 ``` --- ## The consequences of survey design | Design issue | Representative? | Design effect (change to SE) | |:-----------------|:--------------------------------------------------------------------------------|:--------------------------------------------| | Clustering | Yes, if PSUs are proportionally drawn. Otherwise, weighting necessary. | Increases with difference between clusters. | | Stratification | Yes, if strata are sampled with same probability. Otherwise weighting necessary.| Decreases with homogeneity within strata | | Weights | Only if weights are applied. | Increases with the variance of the weights. | --- ## Correcting for survey design -- ### Basic Weighting R has syntax in many commands to apply sampling weights to get representative statistics (e.g. `weighted.mean`, the `weight` option in the `lm` command), but this approach will not correctly adjust standard errors for design effects. -- ### Robust standard errors Robust standard errors will correct standard errors for differential weights, but not for clustering and stratification design effects. -- ### Survey package The `survey` package in R will allow you to specify design and correctly adjust standard errors. --- ## Add Health survey design .pull-left[ ![](images/addhealthdesign.jpeg) ] .pull-right[ - Schools were the primary PSU in a cluster sampling technique, but were also stratified by region, urbanicity, school type, ethnic mix, and size. - Students within schools were stratified by grade and sex and then sampled. - Several oversamples were conducted of ethnic groups and genetically related pairs of students, as well as saturated samples from 16 schools. - Post-stratification adjustments were made to sampling weights to account for region of the country. - The Add Health documentation indicates that the REGION variable should be applied as a stratification variable but this variable is not available in the public release data so we will focus on the design effects of weights and clustering. ] --- ## Sample design variables .pull-left[ ```r head(popularity[,c("cluster","sweight")]) ``` ``` ## cluster sweight ## 1 472 3663.6562 ## 2 472 4471.7647 ## 3 142 3325.3294 ## 4<PHONE_NUMBER> ## 5 145 3939.9402 ## 6<PHONE_NUMBER> ``` - `cluster` is a numeric id indicating the school the student was sampled from. - `sweight` is the inverse of the probability of being sampled. ] .pull-right[  ] --- ## Add Health example, naive approach Lets start with a model that makes no adjustment. ```r model.basic <- lm(indegree~nsports, data=popularity) summary(model.basic)$coef ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3.9969128 0.07221383 55.34830 0.000000e+00 ## nsports 0.5044705 0.04282271 11.78044 1.462129e-31 ``` -- Because I am not using weights, this value will not be representative of all American high school students in 1994-95. I can add the weights provided by Add Health to make this representative: -- ```r model.weight <- update(model.basic, weight=sweight) summary(model.weight)$coef ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.0268786 0.07321643 54.99966 0.00000e+00 ## nsports 0.5330072 0.04347144 12.26109 5.23695e-34 ``` --- ## Add Health example, robust SE approach Adding weights to the `lm` command gives me representative values for the slope and intercept, but my SEs do not correct for the substantial variation in weights. I can correct for this design issue by using robust standard errors: -- ```r library(sandwich) library(lmtest) model.robust <- coeftest(model.weight, vcov=vcovHC(model.weight, "HC1")) model.robust ``` ``` ## ## t test of coefficients: ## ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.026879 0.080106 50.2692 < 2.2e-16 *** ## nsports 0.533007 0.063041 8.4549 < 2.2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` Estimates are identical but standard errors are slightly larger. However, this technique still does not correct for clustering. --- ## Add Health example, using `survey` library I can use the `svydesign` command in the `survey` package to correctly specify both the weights and the clustering in Add Health. The `ids` argument expects a variable name that identifies the clusters by id (in this case, the school id) and the `weight` argument expects a variable name for the weights used. ```r library(survey) popularity.svy <- svydesign(ids=~cluster, weight=~sweight, data=popularity) model.svy <- svyglm(indegree~nsports, design=popularity.svy) summary(model.svy)$coef ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.0268786 0.11483268 35.067357 6.996643e-62 ## nsports 0.5330072 0.07318866 7.282648 5.034252e-11 ``` Estimates are identical to the weighted regression model before, but standard errors have increased due to variation in weights and the cluster design effect. --- ## Comparison of methods .stargazer[

Linear models predicting number of friend nominations by sports played

 	 unweighted 	 weighted 	 robust SE 	 survey: weights 	 survey: weights+cluster

 Intercept	 3.997***	 4.027***	 4.027***	 4.027***	 4.027***

 	 (0.072)	 (0.073)	 (0.080)	 (0.080)	 (0.115)

 Number of sports	 0.504***	 0.533***	 0.533***	 0.533***	 0.533***

 	 (0.043)	 (0.043)	 (0.063)	 (0.063)	 (0.073)

 R2	 0.031	 0.033

 Num. obs.	 4397	 4397	 	 4397	 4397

 ***p < 0.001; **p < 0.01; *p < 0.05

 ] -- - All models except the unweighted version produce the same estimates of slope and intercept based on the weights. -- - The survey weighted and robust SE models both produce the same standard errors. This is because both models account for the design effect of weight heterogeneity but not clustering. --- class: inverse, center, middle background-image: url(images/kyaw-tun-UvL9JrLkagM-unsplash.jpg) background-size: cover name: missing-data # Missing Data --- ## The reality of missing data Missing data exists in most real-world data sets. Therefore, its important to know how to handle missing data in order to know how to properly conduct an analysis. -- .pull-left[ ### Valid Skip ![](images/baby_boss.jpg) Valid skips most arise when a follow-up question is only asked of respondents who gave a certain response to the initial question. If you construct a variable correctly, valid skips should not be considered missing values. ] -- .pull-right[ ### Item non-response ![](images/grumpy_old_man.jpg) Item non-response occurs when respondents fail to respond to a specific question. This may be because they don't know the correct response or they do not feel comfortable answering the question. ] --- ## Example of a valid skip The GSS uses three variables to determine respondents' religious affiliation. - `relig` asks for major religious affiliations such as Catholic, Protestant, Jewish, Muslim, etc. - **If and only if** respondents indicate they are Protestant, they are asked a follow up question recorded in `denom` which asks for their specific denomination. - `denom` only lists major Protestant denominations. If the respondent checks "other", their specific write-in response is recorded in a third variable titled `other`. -- ```r summary(relig[,c("relig","denom","other")]) ``` ``` ## relig denom other ## protestant:34597 other : 7491 pentecostal : 878 ## catholic :14532 baptist-dk which: 6078 mormon : 714 ## none : 6635 southern baptist: 3613 church of christ : 648 ## jewish : 1195 no denomination : 3101 christian : 538 ## other : 1025 united methodist: 2693 jehovah's witnesses: 397 ## (Other) : 1592 (Other) :11959 (Other) : 4224 ## NA's : 23 NA's :24664 NA's :52200 ``` -- There are a lot of missing values for `denom` and `other`, but these are all valid skips based on prior responses. The only true missing values in this set of variables are the 23 respondents who did not respond to the initial question on `relig`. --- ## Kinds of missingness .left-column[ ![](images/coffee_spill.jpg) ] .right-column[ ### Missing Completely at Random (MCAR) Every observation has the same probability of missingness and the missingness of a variable has no relationship to other observed or unobserved variables. If this is true, then removing observations with missing values will not bias results. ### Missing at Random (MAR) The different probabilities of missingness can be fully accounted for by other observed variables in the dataset. If this is true, then various techniques can be used to produce unbiased results by **imputing** values for the missing values. ### Not Missing at Random (NMAR) The different probabilities of missingness depend both on observed and unobserved variables. In this case, we cannot fully correct for bias that might result from missing data. ] ??? In practice, its impossible to distinguish perfectly between MAR and NMAR. If we use many predictors to impute a moderate number of missing values for a case then the MAR assumption is reasonable in the sense that the remaining bias will likely be minimal. --- ## Add Health income example As an example, I will use parental income from the Add Health data to predict popularity. Income is recorded in thousands of dollars, and I have top-coded the values to $200,000. Income is notorious as a variable that will typically have a high non-response rate. The Add Health data are no different: ```r summary(addhealth$parentinc) ``` ``` ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ## 0.00 23.00 40.00 45.37 60.00 200.00 1027 ``` ```r mean(is.na(addhealth$parentinc)) ``` ``` ## [1] 0.2335683 ``` Income is missing for 1027 cases which is roughly a quarter of the dataset. * Note that the `popularity` dataset will not work here as I already imputed missing values. Use the `addhealth` dataset instead. --- ## Two basic approaches to missing values -- .pull-left[ ![](images/throw_away.png) ### Remove cases - Case deletion is easy to implement. - Case deletion assumes MCAR. - Case deletion may result in a substantial reduction is sample size. ] -- .pull-right[ .center[![](images/rabbit.jpg)] ###Impute cases - Imputation is more difficult to implement. - If done properly, imputation assumes MAR. - Imputation will not reduce sample size. ] --- ## Two kinds of case deletion -- .pull-left[ ### Complete-case analysis (listwise deletion) - All cases that have missing values on any of the variables that will be used at some point in the analysis are removed from the start. - This ensures that all results are based on the same sample. ] -- .pull-right[ ### Available-case analysis (pairwise deletion) - Cases are removed model by model or statistic by statistic when the variables used in that particular statistic or model have missing values. - This is by default what will happen in R across different `lm` models with different variables. - This approach allows the researcher to use more data, but different statistics and models will use different subsets of the full data which makes comparability problematic. ] --- ## Removing cases, Add Health example ### Available-case analysis This is what happens by default when you just run nested models in R ```r model1.avail <- lm(indegree~nsports, data=addhealth) model2.avail <- update(model1.avail, .~.+alcoholuse+smoker) model3.avail <- update(model2.avail, .~.+parentinc) ``` -- ### Complete-case analysis To do this, we need to use `na.omit` on the Add Health variables that will be in the most complex model to make sure all models work with the same subset. ```r addhealth.complete <- na.omit(addhealth[,c("indegree","nsports","alcoholuse","smoker","parentinc")]) model1.complete <- lm(indegree~nsports, data=addhealth.complete) model2.complete <- update(model1.complete, .~.+alcoholuse+smoker) model3.complete <- update(model2.complete, .~.+parentinc) ``` --- ## What is different? .pull-left[ .stargazer[ ### Available-case method

Statistical models

 	 Model 1	 Model 2	 Model 3

 Intercept	 4.00***	 3.85***	 3.39***

 	 (0.07)	 (0.08)	 (0.12)

 Number of sports	 0.50***	 0.51***	 0.45***

 	 (0.04)	 (0.04)	 (0.05)

 Drinker	 	 0.71***	 0.66***

 	 	 (0.16)	 (0.18)

 Smoker	 	 0.19	 0.38*

 	 	 (0.16)	 (0.19)

 Parental income (1000s USD)	 	 	 0.01***

 	 	 	 (0.00)

 R2	 0.03	 0.04	 0.05

 Num. obs.	 4397	 4343	 3332

 ***p < 0.001; **p < 0.01; *p < 0.05

 ] ] .pull-right[ .stargazer[ ### Complete-case method

Statistical models

 	 Model 1	 Model 2	 Model 3

 Intercept	 4.05***	 3.88***	 3.39***

 	 (0.08)	 (0.09)	 (0.12)

 Number of sports	 0.49***	 0.49***	 0.45***

 	 (0.05)	 (0.05)	 (0.05)

 Drinker	 	 0.72***	 0.66***

 	 	 (0.18)	 (0.18)

 Smoker	 	 0.37*	 0.38*

 	 	 (0.19)	 (0.19)

 Parental income (1000s USD)	 	 	 0.01***

 	 	 	 (0.00)

 R2	 0.03	 0.04	 0.05

 Num. obs.	 3332	 3332	 3332

 ***p < 0.001; **p < 0.01; *p < 0.05

 ] ] -- ??? - Number of observations changes across models in available-case analysis. - The results for smoking in the complete-case analysis suggest that the difference from model 2 to model 3 in available-case analysis was driven by sample size change, not controlling for income. --- ## Imputation -- .pull-left[ ### Predictive or non-Predictive? - Imputations vary by whether or not they use other predictors in the data to assign imputed values to missing values. - Predictive imputation is done via a statistical model. - Predictive imputation moves the assumption from MAR to MCAR if the model is good. - Non-predictive imputation will systematically bias correlation between variables downward. ] -- .pull-right[ ### Randomness? - Imputations vary by whether or not they include a random component or are purely deterministic. - Deterministic imputation techniques will underestimate variance in the imputed variable and will therefore underestimate standard errors. - The use of randomization leads to another source of variation called **imputation variation** which can only be fully addressed through the technique of **multiple imputation**. ] --- ## Non-Predictive Imputation A very simple (and poor) technique would be just to substitute the mean for valid responses for all missing values. ```r addhealth$parentinc.meani <- ifelse(is.na(addhealth$parentinc), mean(addhealth$parentinc,na.rm=TRUE), addhealth$parentinc) ``` -- Another similar technique that allows for more randomness is just to sample a random valid response on the same variable for each missing value. ```r addhealth$parentinc.randi <- addhealth$parentinc addhealth$parentinc.randi[is.na(addhealth$parentinc)] <- sample(na.omit(addhealth$parentinc), sum(is.na(addhealth$parentinc))) ``` --- ## Mean and random imputation, Add Health .pull-left[  ] .pull-right[  ] --- ## Non-predictive imputation == Bad | Sample | `\(r\)` (indegree and income) | SD (income) | |:--------------------|--------------------------:|-------------------:| | Valid cases | 0.1320575 | 33.6984995| | Valid cases +mean imputed | 0.1167858 | 29.5006936| | Valid cases +random imputed | 0.1076341 | 33.8174488| - Both techniques will systematically underestimate correlation. - Mean imputation will underestimate the variance of the imputed variable. --- ## Quick and dirty method If I do a mean imputation, I can also add the boolean variable indicating missingness as a predictor in my model. When I do this, the effect of the imputed variable will be the same as if if I had thrown out missing values (because we are controlling for missingness), but I can use the full data. ```r summary(lm(indegree~parentinc, data=addhealth))$coef ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3.93119831 0.106835499 36.796742 1.886044e-249 ## parentinc 0.01461554 0.001890363 7.731602 1.392430e-14 ``` ```r summary(lm(indegree~parentinc.meani+is.na(parentinc), data=addhealth))$coef ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3.93119831 0.105945249 37.105942 2.046037e-262 ## parentinc.meani 0.01461554 0.001874611 7.796570 7.888129e-15 ## is.na(parentinc)TRUE -0.18735131 0.130692401 -1.433529 1.517780e-01 ``` - This model still assumes MCAR and reduces variance in the independent variable. - Its primary advantage is that it is a quick method to avoid having to throw out cases that have valid data on other important variables. --- ## Mean imputation with dummy, Add Health  --- ## Regression imputation - I will predict the value of parental income by other independent variables (but never the dependent variable) using a regression model. - In this case, I will transform parental income by the square root as well since it is heavily right skewed. ```r model <- lm(sqrt(parentinc)~race+pseudoGPA+honorsociety+alcoholuse+smoker +bandchoir+academicclub+nsports, data=addhealth) ``` -- Then, I use the `predict` command to get predicted values for all observations and impute the predicted values (or their square, technically) for missing values. ```r predicted <- predict(model, addhealth) addhealth$parentinc.regi <- addhealth$parentinc incmiss <- is.na(addhealth$parentinc) addhealth$parentinc.regi[incmiss] <- predicted[incmiss]^2 summary(addhealth$parentinc.regi) ``` ``` ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ## 0.00 25.19 39.83 43.87 54.20 200.00 39 ``` I still have some missing values, because there were missing values on the variables I used to predict parental income. --- ## Random regression imputation The previous model is deterministic, and will underestimate variance in parental income but I can add a random component to this by sampling from a normal distribution with a mean of zero and standard deviation equal to that of the model residuals. ```r addhealth$parentinc.rregi <- addhealth$parentinc addhealth$parentinc.rregi[incmiss] <- (predicted[incmiss]+ * rnorm(sum(incmiss), 0, sigma(model)))^2 sd(addhealth$parentinc, na.rm=TRUE) ``` ``` ## [1] 33.6985 ``` ```r sd(addhealth$parentinc.regi, na.rm=TRUE) ``` ``` ## [1] 30.10024 ``` ```r sd(addhealth$parentinc.rregi, na.rm=TRUE) ``` ``` ## [1] 32.80286 ``` --- ## Regression imputation, Add Health .pull-left[  ] .pull-right[  ] --- ## Chained equations As you saw above, predictive imputation still leads to some missing values if the variables used for prediction also contain missing values. We can get around this via an iterative procedure called **chained equations**: -- 1. All missing values are given some placeholder value. This might be the mean value, for example. -- 2. For one variable, the placeholder values are removed and missing values put back in. These missing values are then predicted and imputed by some model for this variable. -- 3. Step 2 is repeated for all variables with missing values. When all variables have been imputed, we have completed one iteration. -- 4. Steps 2-3 are then repeated again for some number of iterations. The number of iterations necessary may vary by data, but five iterations is typical. --- ## Using `mice` to impute a single dataset The `mice` command will use chained equations to impute missing values on *all* variables when a dataset is fed in. I can then use the `complete` command to extract a full dataset with no missing values. ```r library(mice) imputed <- mice(addhealth[,c("indegree","race","sex","grade","pseudoGPA","honorsociety", "alcoholuse","smoker","bandchoir","academicclub","nsports", "parentinc")], m=1, print=FALSE) addhealth.ce <- complete(imputed, 1) apply(is.na(addhealth.ce), 2, sum) ``` ``` ## indegree race sex grade pseudoGPA honorsociety ## 0 0 0 0 0 0 ## alcoholuse smoker bandchoir academicclub nsports parentinc ## 0 0 0 0 0 0 ``` I can then use this new dataset `addhealth.ce` to run my models. --- ## Comparison of methods .stargazer[

Models predicting number of friend nominations, Add Health

 	 deletion	 mean	 mean + dummy	 random	 regression	 random regression	 chained equations

 Intercept	 3.391***	 3.360***	 3.399***	 3.470***	 3.350***	 3.441***	 3.340***

 	 (0.120)	 (0.111)	 (0.115)	 (0.104)	 (0.107)	 (0.105)	 (0.103)

 Parental income (1000s)	 0.012***	 0.012***	 0.012***	 0.009***	 0.013***	 0.010***	 0.013***

 	 (0.002)	 (0.002)	 (0.002)	 (0.002)	 (0.002)	 (0.002)	 (0.002)

 Smoker	 0.383*	 0.201	 0.203	 0.196	 0.187	 0.193	 0.201

 	 (0.188)	 (0.161)	 (0.161)	 (0.161)	 (0.162)	 (0.162)	 (0.159)

 Drinker	 0.663***	 0.666***	 0.668***	 0.679***	 0.665***	 0.665***	 0.656***

 	 (0.180)	 (0.155)	 (0.155)	 (0.155)	 (0.155)	 (0.156)	 (0.154)

 Number of sports	 0.454***	 0.476***	 0.474***	 0.481***	 0.462***	 0.475***	 0.461***

 	 (0.049)	 (0.043)	 (0.043)	 (0.043)	 (0.043)	 (0.043)	 (0.043)

 Income missing	 	 	 -0.162

 	 	 	 (0.130)

 R2	 0.048	 0.046	 0.046	 0.044	 0.048	 0.045	 0.049

 Num. obs.	 3332	 4343	 4343	 4343	 4320	 4320	 4397

 ***p < 0.001; **p < 0.01; *p < 0.05

 ] --- ## Imputation variability -- - The methods of random imputation have the benefit of preserving the standard deviation of the imputed variable and therefore calculating correct standard errors, but they also introduce a new source of uncertainty. -- - Each time I do an imputation with a random component (e.g. random regression, chained equation), I will get a somewhat different set of values. -- - Therefore, we now have **imputation variability** to add to our inferential concerns alongside **sampling variability**. -- - We can use **multiple imputation** to adjust our results for **imputation variability**. --- ## Multiple imputation process -- 1. Use imputation process with random component to impute missing values and repeat this process to produce `\(m\)` separate **complete datasets**. Each of these datasets will be somewhat different due to the randomization of imputation. Usually `\(m=5\)` is sufficient. -- 2. Run `\(m\)` separate parallel models on each imputed **complete dataset**. As a result, you will have `\(m\)` sets of regression coefficients and standard errors. -- 3. Pool the regression coefficients across datasets by taking the mean across all `\(m\)` datasets. -- 4. Pool standard errors by taking the mean across all `\(m\)` datasets *plus* the between model standard deviation in coefficients. The formula for the overall standard error is: `$$SE_{\beta}=\sqrt{W+(B+\frac{B}{m})}$$` Where `\(W\)` is the squared mean standard error across all `\(m\)` datasets, and `\(B\)` is the variance in coefficient estimates calculated across all `\(m\)` models. --- ## Multiple imputation with the `mice` package To get multiple complete datasets with `mice` we just need to specify the number of distinct complete datasets to create in the second argument: ```r imputations <- mice(popularity, 5, printFlag=FALSE) ``` -- - The imputations object now contains five fully complete imputed datasets. Its possible to extract any one of these datasets with the command `complete(imputations, i)` where `i` is replaced by a number between 1 and 5. -- - We can now conduct our parallel analysis on these five datasets and combine results. There is easy and a hard way to do this. Its useful to know both. --- ## Easy way: let `mice` do the hard work The `mice` package has a lot of nice features, including an object specific function for the `with` command and a `pool` command that make multiple imputation as easy as falling off a log: ```r model.mi <- pool(with(imputations, lm(indegree~parentinc+smoker+alcoholuse+nsports))) summary(model.mi) ``` ``` ## estimate std.error statistic df p.value ## (Intercept) 3.38654549 0.102909528 32.907988 4389.514 0.000000e+00 ## parentinc 0.01126339 0.001630526 6.907827 4389.514 5.627498e-12 ## smokerSmoker 0.23538273 0.159687599 1.474020 4389.514 1.405479e-01 ## alcoholuseDrinker 0.64050221 0.154054327 4.157639 4389.514 3.276895e-05 ## nsports 0.47096906 0.042956107 10.963961 4389.514 0.000000e+00 ``` --- ## Hard way: for-loop The hard way isn't really that hard. It is useful to know for cases where the easy way won't work, such as when you need to run complicated models (using `svydesign` would be an example). ```r b <- se <- NULL for(i in 1:5) { imputation <- complete(imputations,i) model <- lm(indegree~parentinc+smoker+alcoholuse+nsports, data=imputation) b <- cbind(b, coef(model)) se <- cbind(se, summary(model)$coef[,2]) } b ``` ``` ## [,1] [,2] [,3] [,4] [,5] ## (Intercept) 3.38654549 3.38654549 3.38654549 3.38654549 3.38654549 ## parentinc 0.01126339 0.01126339 0.01126339 0.01126339 0.01126339 ## smokerSmoker 0.23538273 0.23538273 0.23538273 0.23538273 0.23538273 ## alcoholuseDrinker 0.64050221 0.64050221 0.64050221 0.64050221 0.64050221 ## nsports 0.47096906 0.47096906 0.47096906 0.47096906 0.47096906 ``` The `b` and `se` objects are matrices that contain the coefficients and standard errors, respectively, for each model on the column. --- ## Hard way: pool the results Now we can pool the results using the `b` and `se` matrices and some creative use of `apply` commands. ```r b.pool <- apply(b,1,mean) between.var <- apply(b,1,var) within.var <- apply(se^2,1,mean) se.pool <- sqrt(within.var+between.var+between.var/5) t.pool <- b.pool/se.pool pvalue.pool <- (1-pnorm(abs(t.pool)))*2 data.frame(b.pool, se.pool, t.pool, pvalue.pool) ``` ``` ## b.pool se.pool t.pool pvalue.pool ## (Intercept) 3.38654549 0.102909528 32.907988 0.000000e+00 ## parentinc 0.01126339 0.001630526 6.907827 4.921397e-12 ## smokerSmoker 0.23538273 0.159687599 1.474020 1.404762e-01 ## alcoholuseDrinker 0.64050221 0.154054327 4.157639 3.215542e-05 ## nsports 0.47096906 0.042956107 10.963961 0.000000e+00 ``` --- class: inverse, center, middle background-image: url(images/nareeta-martin-z7LRp54Yulw-unsplash.jpg) background-size: cover name: multicollinearity # Multicollinearity and Scale Creation --- ## The problem of multicollinearity - Given that adding more independent variables to your model allows you to account for potential omitted variable bias, why wouldn't you just put in as many variables as you can? -- - Because of **multicollinearity**. Multicollinearity occurs when there is collectively a high correlation between the independent variables in the model. -- - Intuitively, its easy to understand why multicollinearity is a problem. When two independent variables are highly correlated with one another, it becomes hard to separate out their effect on the dependent variable. -- - Technically, the effect of multicollinearity is to inflate standard errors and make coefficient estimates highly unstable across different combinations of highly collinear terms in the model. --- ## Two kinds of multicollinearity -- .pull-left[ ### Structural multicollinearity Structural multicollinearity occurs when one independent variable is completely determined by another independent variable or set of independent variables. This really isn't an issue with the data, but rather a specification error by the researcher. ] -- .pull-right[ ### Data-based multicollinearity Data-based multicollinearity occurs when a set of variables are highly but not perfectly correlated with one another in the empirical data. It is a hazard of using observational data where the characteristics for which you want to get separate effects are difficult to separate. ] --- ## Example of structural multicollinearity Lets predict violent crime by % female and % male. ```r crimes$PctFemale <- 100-crimes$PctMale summary(lm(Violent~PctMale+PctFemale, data=crimes)) ``` ``` ## ## Call: ## lm(formula = Violent ~ PctMale + PctFemale, data = crimes) ## ## Residuals: ## Min 1Q Median 3Q Max ## -293.57 -104.01 -54.44 88.98 783.97 ## *## Coefficients: (1 not defined because of singularities) ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4317.26 1628.20 2.652 0.0108 * ## PctMale -79.75 33.01 -2.416 0.0195 * *## PctFemale NA NA NA NA ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 186.5 on 49 degrees of freedom ## Multiple R-squared: 0.1064, Adjusted R-squared: 0.08819 ## F-statistic: 5.836 on 1 and 49 DF, p-value: 0.01948 ``` --- ## Singularity! Not as cool as it sounds. One of the terms was dropped from the model because the terms are perfectly collinear. ```r cor(crimes$PctFemale,crimes$PctMale) ``` ``` ## [1] -1 ``` -- This is not a problem of the model, but our thinking. Either term by itself will give you full information. ```r coef(lm(Violent~PctMale, data=crimes)) ``` ``` ## (Intercept) PctMale ## 4317.26466 -79.75285 ``` ```r coef(lm(Violent~PctFemale, data=crimes)) ``` ``` ## (Intercept) PctFemale ## -3658.02067 79.75285 ``` --- ## Detecting data-based multicollinearity .left-column[ ![](images/detective.jpg) ] .right-column[ - Standard errors increase substantially across nested models with more covariates - Regression coefficients are highly unstable across nested models - Examination of correlation matrix - Calculate variance inflation factors ] --- ## Multicollinearity example: NYC non-profits .pull-left[ We will look at data collected by myself and Nicole Marwell on the spatial distribution of money contracted out to organization for social services by the City of New York from 1997-2001. The unit of analysis is a NYC health area, which can loosely be thought of as a neighborhood. The variables are: - **amountcapita**: The dollar amount of money provided to the health area divided by the population size of the health area. We will log it due to skewness. - **poverty**: percent of population below the poverty line. - **unemployed**: unemployment rate. - **income**: median household income. We will log it due to skewness. ] .pull-right[  ] --- ## Model results, NYC example .right-column[ .stargazer[

Models predicting social service funding per capita (logged)

 	 Model 1	 Model 2	 Model 3	 Model 4	 Model 5	 Model 6	 Model 7

 intercept	 4.14***	 4.13***	 10.24***	 4.12***	 -24.70***	 3.08	 -29.60***

 	 (0.17)	 (0.23)	 (2.05)	 (0.23)	 (4.94)	 (3.27)	 (5.31)

 poverty rate	 0.04***	 	 	 0.04***	 0.13***	 	 0.12***

 	 (0.01)	 	 	 (0.01)	 (0.02)	 	 (0.02)

 unemployment rate	 	 0.09***	 	 0.00	 	 0.10**	 0.08*

 	 	 (0.02)	 	 (0.03)	 	 (0.04)	 (0.03)

 median income (logged)	 	 	 -0.49*	 	 2.52***	 0.09	 2.91***

 	 	 	 (0.19)	 	 (0.43)	 (0.28)	 (0.46)

 R2	 0.08	 0.04	 0.02	 0.08	 0.16	 0.04	 0.18

 Num. obs.	 341	 341	 341	 341	 341	 341	 341

 ***p < 0.001; **p < 0.01; *p < 0.05

 ] ] .left-column[ The models exhibit classic signs of multicollinearity: - Highly unstable coefficients across models. - Large increase in standard errors. ] --- ## The correlation matrix as diagnostic .pull-left[ ```r cor(nyc[,c("poverty","unemployed","lincome")]) ``` ``` ## poverty unemployed lincome ## poverty 1.0000000 0.6982579 -0.9120680 ## unemployed 0.6982579 1.0000000 -0.7419472 ## lincome -0.9120680 -0.7419472 1.0000000 ``` The correlation between the three variables is very high, suggesting strong multicollinearity. Note that, while the correlation matrix is often helpful, it may not reveal the full extent of multicollinearity because it only looks at bivariate relationships between the variables. ] -- .pull-right[ ### Visualize as correlogram ```r library(corrgram) corrgram(nyc[,c("poverty","unemployed", "lincome")], upper.panel="panel.cor", lower.panel="panel.pts") ```  ] --- ## Variance inflation factors The **variance inflation factor** (VIF) is the multiplicative factor by which the variance in the estimation of any one coefficient from the current model is increased due to multicollinearity relative to a model in which the variables are uncorrelated. The square root of the VIF is roughly the expected factor increase in the standard error. It can be shown that the VIF for the `\(i\)`th variable is given by: `$$VIF_i=\frac{1}{1-R_i^2}$$` Where `\(R_i^2\)` is the r-squared value when the given independent variable is predicted by all of the other independent variables in the model. For example, we could calculate the VIF for poverty in the full model by: ```r 1/(1-summary(lm(poverty~unemployed+lincome, data=nyc))$r.squared) ``` ``` ## [1] 5.984485 ``` The square root of this VIF is 2.45, indicating that the standard error for poverty is more than doubled due to multicollinearity in the full model. --- ## Estimating VIF with the `vif` command We can also use the `vif` function in the `car` package to quickly estimate all VIFs for a given model: ```r library(car) vif(model.allthree) ``` ``` ## poverty unemployed lincome ## 5.984485 2.238379 6.822173 ``` The general rule of thumb is that a VIF over four is problematic enough that it needs to be addressed in some manner. Clearly, multicollinearity is a big problem here. --- ## You have multicollinearity, now what? -- ### Remove variables A simple approach is to remove some of the highly correlated covariates. However, because the variables are not perfectly correlated, you are basically throwing out some information. -- ### Separate models Another approach is to run separate models with only one of the highly collinear variables in each model. This can also be unsatisfying because each of the models is underestimating the total effect of the variables collectively. -- ### Create a scale In cases where the collinear variables are all thought to represent the same underlying conceptual variable, we can combine them into a single synthetic scale. --- ## Standardization and reverse ordering .pull-left[ - In many cases, the variables that are thought to make up a scale might be measured in the same manner. - When variables are measured differently (as in our example), then they must be standardized in some way to make them comparable before scale construction and evaluation. The most common way to do this is by creating **z-scores** by subtracting by the mean and dividing by the standard deviation of each variable. - Some variables may be positively related to the underlying concept while others may be negatively related. It may be necessary to reverse the coding of a variable to make all variables positively related to the underlying concept. ] .pull-right[ ### NYC example ```r nyc$poverty.z<-(nyc$poverty-mean(nyc$poverty))/ sd(nyc$poverty) nyc$unemployed.z<-(nyc$unemployed-mean(nyc$unemployed))/ sd(nyc$unemployed) nyc$lincome.z<- -1*(nyc$lincome-mean(nyc$lincome))/ sd(nyc$lincome) ``` Note that I am multiplying the `lincome` result by -1 to reverse code it. I can also use the `scale` command to do this: ```r nyc$poverty.z <- scale(nyc$poverty) nyc$unemployed.z <- scale(nyc$unemployed) nyc$lincome.z <- -1*scale(nyc$lincome) ``` ] --- ## Cronbach's alpha Cronbach's `\(\alpha\)` is a statistic developed in psychology to test the degree to which different variables measure the same underlying concept. -- - Cronbach's `\(\alpha\)` is often thought of as a test of the **internal reliability** of a scale. -- - You can think about Cronbach's `\(\alpha\)` as a summary measure of the correlation matrix we saw earlier. It goes from 0 to 1, with 0 indicating no shared correlation, and 1 indicating perfect correlation between all items. -- The `pysch` package includes and `alpha` command that will calculate `\(\alpha\)`. ```r psych::alpha(nyc[,c("poverty.z","unemployed.z","lincome.z")])$total ``` ``` ## raw_alpha std.alpha G6(smc) average_r S/N ase mean ## 0.9159293 0.9159293 0.9013009 0.784091 10.89474 0.008313234 -2.211598e-16 ## sd median_r ## 0.9252355 0.7419472 ``` The `\(\alpha\)` of 0.915 indicates a very high level of shared covariance between the three variables. --- ## Summated scale Once I have standardized and (if necessary) reverse coded all of the variables for my scale, I can create a simple **summated scale** by adding them up. I am also going to scale this variable with the `scale` command so that it has a mean of zero and a standard deviation of one. ```r nyc$deprivation_summated <- scale(nyc$poverty.z+nyc$unemployed.z+nyc$lincome.z) model <- lm(log(amtcapita)~deprivation_summated, data=nyc) summary(model)$coef ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.9474224 0.08807143 56.175114 8.329185e-174 ## deprivation_summated 0.3747882 0.08820085 4.249258 2.775153e-05 ``` ```r exp(coef(model)) ``` ``` ## (Intercept) deprivation_summated ## 140.811538 1.454683 ``` A one standard deviation increase on my deprivation scale is associated with a 45% increase in the amount for social services in a neighborhood. --- ## Factor analysis Another approach to creating a scale is **factor analysis**. Factor analysis is a method to extract the underlying latent variables (the factors) from a set of observed variables. We will start with an example where we assume a single underlying latent factor. -- Lets say we have three variables `\(z_1\)`, `\(z_2\)`, and `\(z_3\)` all measured as z-scores. We can construct a system of equations that relate each of these variables to a common shared factor `\(F\)` (also on a standard scale with a mean of 0 and SD of 1) and three unique factors, `\(Q_1\)`, `\(Q_2\)`, and `\(Q_3\)`. `$$z_{i1}=b_1F_i+u_1Q_{i1}$$` `$$z_{i2}=b_2F_i+u_2Q_{i2}$$` `$$z_{i3}=b_3F_i+u_3Q_{i3}$$` The `\(F_i\)` are called the **factor scores** and are the values of the common factor for each observation. The `\(b\)` values are called the **factor loadings** and give the correlation between each observed variable and the common factor. The unique factor components are not actually estimated in factor analysis but serve as the "residual" components. --- ## Estimation in factor analysis .pull-left[ The factor scores and factor loadings for a factor analysis can be estimated by mathematical manipulation of the observed correlation matrix between variables, but we won't get into the mathematical details here. - The `factanal` command in base R will estimate a factor analysis by maximum likelihood estimation (a technique we will learn in the next module). - The `fa` command in the `psych` package will also estimate factor analysis using a variety of techniques. It also has some nice additional tools that will make it more useful. ] -- .pull-right[ ### NYC example ```r factor_nyc <- fa(nyc[,c("poverty", "unemployed", "lincome")], 1, rotate="oblimin") ``` - The second argument specifies that we want only one common factor. - We also need to specify a technique of *rotation* for the factor loadings because there are an infinite number of possible ways to express them. The oblimin method is a standard approach that helps to maximize differences between factors without forcing them to be uncorrelated. ] --- ## Factor loadings .pull-left[ ### Extract loadings ```r loadings(factor_nyc) ``` ``` ## ## Loadings: ## MR1 ## poverty 0.926 ## unemployed 0.754 ## lincome -0.984 ## ## MR1 ## SS loadings 2.396 ## Proportion Var 0.799 ``` The bottom part shows that 79.9% of the variation in the three variables is accounted for by the common factor. The loadings themselves show the correlation coefficient between each observed variable and the common factor. ] -- .pull-right[ ### Visualize factor loadings ```r fa.diagram(factor_nyc) ```  ] --- ## Factor Scores I can easily extract my factor scores from the the factor analysis object and use them as my scale measure of deprivation. ```r nyc$deprivation_factor <- factor_nyc$scores[,"MR1"] model <- lm(log(amtcapita)~deprivation_factor, data=nyc) summary(model)$coef ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.9474224 0.08908910 55.533419 2.797390e-172 ## deprivation_factor 0.2849722 0.09036319 3.153631 1.756654e-03 ``` ```r exp(coef(model)) ``` ``` ## (Intercept) deprivation_factor ## 140.811538 1.329725 ``` A one standard deviation increase on my deprivation scale is associated with a 33% increase in the amount for social services in a neighborhood. --- ## Whats the difference? .pull-left[  ] -- .pull-right[ - The summated scale uses all the variation in the three variables to generate the score. - The factor analysis only uses the shared component to generate the score. The unique component is left out. - A closely related technique to factor analysis called **principal component analysis** uses all of the variation in the items and will produce results virtually identical to the summated score. ] --- ## Factor analysis with multiple factors The previous example only used one factor and produced results that were very similar to a simple summated scale. However, it is also possible to use factor analysis to identify more than one common factor shared among a set of variables. -- The formulas for factor analysis above can be generalized to a set of `\(J\)` observed variables and `\(m\)` factors by a set of `\(J\)` equations. For the `\(j\)`th observed variable: `$$z_{ji}=b_{j1}F_{1i}+b_{j2}F_{2i}+\ldots+b_{jm}F_{ji}+u_jQ_{ji}$$` There will be a set of `\(J\)` factor loadings for each of the `\(m\)` factors. The key question with this technique is what are the appropriate number of factors? This is typically determined by an analysis of: -- - how much total variation is explained by a given number of factors. -- - whether the factor loadings for a given number of factors make theoretical sense. --- ## Example: social conservatism among Muslims .pull-left[ Between 2008-12, The Pew Research Center surveyed Muslims from numerous countries around the world and (in addition to other questions) asked respondents about the moral acceptability of these practices: - divorce - polygamy - fertility control - alcohol - euthanasia - suicide - abortion - prostitution - premarital sex - homosexuality ] -- .pull-right[ The data are coded as an ordinal response: - Morally acceptable - Depends/Not a moral issue - Morally wrong I re-code the data on a numeric scale from 1 to 3. ```r morality <- cbind(divorce=pew$moral_divorce, fertility=pew$moral_fertility, alcohol=pew$moral_alcohol, euthansia=pew$moral_euthansia, suicide=pew$moral_suicide, abortion=pew$moral_abortion, prostitution=pew$moral_prostitution, premar_sex=pew$moral_premar_sex, moral_gay=pew$moral_gay) ``` ] --- ## Correlogram of responses ```r corrgram(morality, upper.panel="panel.cor", order="PCA") ```  --- ## Using the correlation matrix as data -- - Factor analysis only needs the correlation matrix between variables rather than the original data itself. -- - This is one case where available-case analysis makes sense. I can maximize the number of cases on each pairwise correlation coefficient to use all of the available data. -- The `use="pairwise.complete.obs"` argument in the `cor` function will use available-case analysis when it computes the correlation coefficient between each pair of variables. ```r morality_r <- cor(morality, use="pairwise.complete.obs") ``` -- ### Try factor analysis with up to three factors ```r morality_fa1 <- fa(morality_r, 1, rotate="oblimin") morality_fa2 <- fa(morality_r, 2, rotate="oblimin") morality_fa3 <- fa(morality_r, 3, rotate="oblimin") ``` --- ## Comparing models with `fa.diagram`  --- ## Extending factor analysis -- Factor analysis is part of a larger family of methods that are often described as **data reduction** techniques: How can you reduce the amount of information in a given set of variables into a smaller set of variables? -- .pull-left[ ### Latent class analysis An important cousin of factor analysis that can be used for categorical variables. ] -- .pull-right[ ### Structural equation modeling The idea of factor analysis can be generalized to construct entire systems of equations that involve both latent and observed variables. ] --- class: inverse, center, middle background-image: url(images/caleb-jones-J3JMyXWQHXU-unsplash.jpg) background-size: cover name: model-selection # Model Selection --- ## Predicting violent crime rates Lets predict violent crime rates in US States by a variety of characteristics. Note that some of these models are **nested** within other models, meaning that the more complex models just add additional variables to earlier models. Other models are non-nested. ```r model.demog <- lm(Violent~MedianAge+PctMale, data=crimes) model.ineq <- update(model.demog,.~.+Gini) model.pov <- update(model.ineq,.~.+Poverty) model.unemp <- update(model.pov,.~.+Unemployment) model.justunemp <- lm(Violent~Unemployment, data=crimes) model.another <- update(model.justunemp,.~.+Gini+MedianAge) ``` --- ## Which model is best? .stargazer[

Models predicting violent crime rate

 	 Model 1	 Model 2	 Model 3	 Model 4	 Model 5	 Model 6

 Intercept	 7534.81***	 619.07	 604.71	 128.07	 -6.09	 -638.35

 	 (1724.03)	 (2995.68)	 (3031.21)	 (2972.28)	 (113.33)	 (599.66)

 Median Age	 -38.56***	 -27.84*	 -28.06*	 -27.84*	 	 -24.90*

 	 (10.89)	 (10.94)	 (11.28)	 (11.02)	 	 (9.40)

 Percent Male	 -115.67***	 -23.28	 -23.21	 -12.29

 	 (31.39)	 (44.75)	 (45.23)	 (44.59)

 Gini	 	 43.01**	 43.71*	 39.04*	 	 38.68**

 	 	 (15.68)	 (17.36)	 (17.15)	 	 (11.67)

 Poverty Rate	 	 	 -0.88	 -6.09

 	 	 	 (8.91)	 (9.17)

 Unemployment Rate	 	 	 	 25.47	 45.89***	 23.02

 	 	 	 	 (14.18)	 (12.99)	 (13.14)

 R2	 0.29	 0.39	 0.39	 0.43	 0.20	 0.42

 Num. obs.	 51	 51	 51	 51	 51	 51

 ***p < 0.001; **p < 0.01; *p < 0.05

 ] --- ## A tradeoff: accuracy vs. parsimony .pull-left[ .center[![](images/bullseye.jpg)] ### Accuracy We want our models to be accurate, or more formally we want them to have high **goodness of fit**. This is measured in OLS regression models primarily by the `\(R^2\)` value. A more accurate model accounts for more of the variation in the dependent variable. ] -- .pull-right[ .center[![](images/piggybank.jpg)] ### Parsimony We want our models to be parsimonious. They should account for variation in the dependent variable with a minimal number of explanatory variables. Models that are not parsimonious are rarely theoretically very interesting. ] ??? We can always add more variables, interactions, non-linear terms, and so forth to get a more accurate model, but this comes at the cost of parsimony. Ultimately, this would lead to us simply throwing up our hands and saying "everybody is different." We knew that already. Our goal is to see if a substantial part of the observed variation in outcome can be explained by a minimal set of predictors, in theoretically expected (or sometimes unexpected) ways. --- ## How to choose a model -- .pull-left[ ### Do Nots - **Do not** choose a model based on all of the coefficients being statistically significant. That is a misuse of statistical inference based on an incorrect interpretation of a p-value. - **Do not** choose a model because it produces the results you like the best. That is not science. - **Do not** choose a model based on some purely mechanical algorithm. Model choice needs to be driven by what questions you want to ask as much as goodness of fit and parsimony. - **Do not** choose a model based on `\(R^2\)` values. An additional variable will make the `\(R^2\)` value somewhat better, so there is no end to the complexity of the model you fit. ] -- .pull-right[ ### Good Practices - **Do** think theoretically about what is the best model for the question at hand. - **Do** make use of model selection tools to help make a decision, but not in a mechanistic way. Use multiple tools. - **Do** conduct and report sensitivity analyses showing how dependent your result is upon the selected model vs. other models. - **Do** weigh both parsimony and accuracy, relative to the goals of your research question. ] ??? Maybe the control variable needs to be in there because prior work or theory says its important, regardless of what the model fit statistics say. Perhaps interactions need to be in the model because your key interest is in intersectionality. And so on. --- ## Model selection tools - The F-Test - Adjusted `\(R^2\)` - Bayesian Information Criterion (BIC) - Bayesian Model Averaging (BMA) --- ## Null vs. saturated model Many model selection tools utilize a distinction between the **null model** and the **saturated model**. -- .pull-left[ ### Null model `$$\hat{y}_i=\bar{y}$$` The null model is the regression model with no predictors. This model is the most parsimonious, but the least accurate. ] -- .pull-right[ ### Saturated model `$$\hat{y}_i=\beta_0+\beta_1x_{i1}+\beta_2x_{i2}+\ldots+\beta_px_{ip}$$` Where, `$$n=p$$` The saturated model is a hypothetical model in which the number of predictors equals the number of observations and we fit each `\(y_i\)` exactly. This model is the most accurate with `\(R^2=1\)`, but the least parsimonious. ] --- ## The F-Test The F-test is a classical hypothesis test for comparing two different models in which one model is nested within the other. This means that the second model includes all of the same independent variables as model 1 as well as an additional `\(g\)` independent variables. The null hypothesis for the F-test is that the effects for all of the additional `\(g\)` variables are zero. -- The F-test produces a test statistic called the **F-statistic**: `$$F=\frac{(SSR_1-SSR_2)/g}{SSR_2/(n-g-k-1)}$$` Where `\(SSR\)` is the sum of squared residuals for a given model, `\(g\)` is the additional terms added to the second model, and `\(k\)` is the number of terms in the first model. The F-statistic is a ratio of the mean amount of variation explained by the new predictors in the second model (numerator) by the mean amount of variation unexplained per degree of freedom in the second model (denominator). -- Assuming the null hypothesis is true, the F-statistic will come from an F-distribution with degrees of freedom equal to `\(g\)` and `\(n-g-k-1\)` (the F-distribution has two DF parameters). Thus, its possible to determine how far out in the tail of this distribution the estimated F-statistic is to calculate a p-value for the hypothesis test. --- ## Calculating F-statistic .pull-left[ ```r model.null <- lm(Violent~1, data=crimes) SSR1 <- sum(model.null$residuals^2) SSR2 <- sum(model.demog$residuals^2) g <- length(model.demog$coef)-length(model.null$coef) k <- length(model.null$coef)-1 n <- length(model.null$residuals) Fstat <- ((SSR1-SSR2)/g)/(SSR2/(n-g-k-1)) Fstat ``` ``` ## [1] 9.869375 ``` ```r 1-pf(Fstat, g, n-g-k-1) ``` ``` ## [1] 0.0002568606 ``` The p-value for the hypothesis test is 0.00026. So, I would likely not reject the null hypothesis and would prefer the null model to the model predicting violent crime by median age and percent male. ] .pull-right[  ] --- ## Using `anova` command for F-test Rather than do this by hand, you can feed both models into the `anova` command. ```r anova(model.null, model.demog) ``` ``` ## Analysis of Variance Table ## ## Model 1: Violent ~ 1 ## Model 2: Violent ~ MedianAge + PctMale ## Res.Df RSS Df Sum of Sq F Pr(>F) ## 1 50 1907432 ## 2 48 1351616 2 555817 9.8694 0.0002569 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- ## Multiple model comparisons The `anova` command is not limited to two models. We can feed in a list of models **so long as they are nested models.** ```r anova(model.null, model.demog, model.ineq, model.unemp) ``` ``` ## Analysis of Variance Table ## ## Model 1: Violent ~ 1 ## Model 2: Violent ~ MedianAge + PctMale ## Model 3: Violent ~ MedianAge + PctMale + Gini ## Model 4: Violent ~ MedianAge + PctMale + Gini + Poverty + Unemployment ## Res.Df RSS Df Sum of Sq F Pr(>F) ## 1 50 1907432 ## 2 48 1351616 2 555817 11.5056 9.209e-05 *** ## 3 47 1165120 1 186496 7.7211 0.007935 ** ## 4 45 1086941 2 78179 1.6183 0.209554 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` The F-statistic is somewhat different than the two model comparison because the denominator for each F-statistic always uses the most complex model. --- ## F-test for adding one variable ```r anova(model.demog, model.ineq) ``` ``` ## Analysis of Variance Table ## ## Model 1: Violent ~ MedianAge + PctMale ## Model 2: Violent ~ MedianAge + PctMale + Gini ## Res.Df RSS Df Sum of Sq F Pr(>F) ## 1 48 1351616 ## 2 47 1165120 1 186496 7.5231 0.008593 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ```r summary(model.ineq)$coef["Gini",] ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## 43.006211344 15.679510943 2.742828619 0.008593105 ``` The p-value for an F-test where only one variable is added is identical to the p-value for a t-statistic of that particular regression coefficient. So, the F-test really only adds something when comparing across models where the more complex model adds multiple new variables. --- ## Limitations of the F-test -- - The F-test can only be used to compare nested models. This limits its use substantially because we can imagine that it might be theoretically important to compare non-nested models. -- - The F-test for adding a single variable is equal to the same t-test for the regression coefficient of that variable in the new model. If we just add each new variable sequentially to our nested models, the F-test adds nothing new. If we add in clusters of variables, the F-test does not help us distinguish which of the variables in the cluster were most important. -- - The F-test is based on the logic of null hypothesis significance testing. That means it is largely driven by sample size. The preferred model by the F-test will differ across two different datasets if they vary in sample size but are the same otherwise. --- ## Methods for penalizing goodness of fit -- - `\(R^2\)` can be thought of as a goodness-of-fit statistic. It measures the accuracy of the model in predicting the dependent variable. -- - `\(R^2\)` cannot be used directly to select models because it can only get bigger as more independent variables are added to the model. It will always prefer the more complex model unless there is literally zero correlation between the new independent variables and the dependent variable. -- - What if we could penalize `\(R^2\)` by some value that depends on the number of independent variables in the model? This would give us a single measure that combines accuracy and parsimony. --- ## Adjusted `\(R^2\)` The formula for adjusted `\(R^2\)` is: `$$R^2_{adj}=R^2-(\frac{p}{n-p-1})(1-R^2)$$` Where `\(p\)` is the number of independent variables in the model. The subtracted part is the parsimony penalty and it is a function of the number of variables as a proportion of the sample size and the lack of goodness-of-fit, `\(1-R^2\)`. -- This value is calculated by default by the `summary` command for a linear model in R. It also shows up in the default output to `texreg`. ```r summary(model.demog)$r.squared ``` ``` ## [1] 0.2913952 ``` ```r summary(model.demog)$adj.r.squared ``` ``` ## [1] 0.26187 ``` --- ## Bayesian Information Criterion There are two common forms of BIC: - BIC is implicitly compared to the saturated model. The `BIC` command in R will give you this value. - BIC' is implicitly compared to the null model. -- The equation for BIC' for OLS regression models is most intuitive: `$$BIC'=n\log(1-R^2)+p\log(n)$$` Where `\(p\)` is the number of independent variables in the model. The first component is the accuracy of the model and the second component is the parsimony penalty. -- BIC is unusual in that lower values are better, and negative values are most preferred. A negative value means you prefer this model to its comparison (e.g. null or saturated model). You can directly compare the values of any two given models by BIC' or BIC, and you will prefer the model with the lower value. The models do not need to be nested. --- ## Calculating BIC ```r bic.null <- function(model) { rsq <- summary(model)$r.squared n <- length(model$resid) p <- length(model$coef)-1 return(n*log(1-rsq)+p*log(n)) } bic.null(model.demog) ``` ``` ## [1] -9.703675 ``` ```r BIC(model.demog) ``` ``` ## [1] 679.8933 ``` ```r BIC(model.demog)-BIC(model.null) ``` ``` ## [1] -9.703675 ``` --- ## BIC or BIC' give the same difference ```r bic.null(model.unemp)-bic.null(model.demog) ``` ``` ## [1] 0.6808653 ``` ```r BIC(model.unemp)-BIC(model.demog) ``` ``` ## [1] 0.6808653 ``` -- We slightly prefer the more complex model with unemployment, poverty, and gini in addition to the demographic characteristics. -- A general rule of thumb is that BIC differences less than two give weak preferences, while BIC differences greater than six give strong preferences, and BIC differences greater than ten give very strong preferences. --- ## Model comparison .stargazer[

Models predicting violent crime rate

 	 Model 1	 Model 2	 Model 3	 Model 4	 Model 5	 Model 6

 Intercept	 7534.81***	 619.07	 604.71	 128.07	 -6.09	 -638.35

 	 (1724.03)	 (2995.68)	 (3031.21)	 (2972.28)	 (113.33)	 (599.66)

 Median Age	 -38.56***	 -27.84*	 -28.06*	 -27.84*	 	 -24.90*

 	 (10.89)	 (10.94)	 (11.28)	 (11.02)	 	 (9.40)

 Percent Male	 -115.67***	 -23.28	 -23.21	 -12.29

 	 (31.39)	 (44.75)	 (45.23)	 (44.59)

 Gini	 	 43.01**	 43.71*	 39.04*	 	 38.68**

 	 	 (15.68)	 (17.36)	 (17.15)	 	 (11.67)

 Poverty Rate	 	 	 -0.88	 -6.09

 	 	 	 (8.91)	 (9.17)

 Unemployment Rate	 	 	 	 25.47	 45.89***	 23.02

 	 	 	 	 (14.18)	 (12.99)	 (13.14)

 R2	 0.29	 0.39	 0.39	 0.43	 0.20	 0.42

 Adjusted R2	 0.26	 0.35	 0.34	 0.37	 0.19	 0.39

 BIC (null)	 -9.70	 -13.34	 -9.42	 -9.02	 -7.63	 -16.28

 ***p < 0.001; **p < 0.01; *p < 0.05

 ] ??? BIC prefers model 6, while adjusted `\(R^2\)` prefers model 7 and is indifferent between models 5 and 6. --- ## Model uncertainty -- - Given a set of potential independent variables, there is a very large number of potential models that could be chosen, even if we don't consider added complications like interaction or polynomial terms. -- - Generally, we select only one model from this full set. Even when we examine alternative models, we don't examine every possible set.  -- - The `BMA` package in R will do **Bayesian Model Averaging** which will provide an expected value for every coefficient across all possible models, as well as the probability a variable is in the model, and the top preferred models. --- ## Bayesian Model Averaging ```r library(BMA) model.bma <- bic.glm(crimes[,c("MedianAge","PctMale","PctLessHS","MedianIncomeHH", "Unemployment","Poverty", "Gini")], crimes$Violent, glm.family=gaussian) summary(model.bma) ``` ``` 29 models were selected Best 5 models (cumulative posterior probability = 0.5217 ): p!=0 EV SD model 1 model 2 model 3 model 4 model 5 Intercept 100 -54.3247 858.601 58.711 9.610 565.958 112.820 432.648 MedianAge 25.8 -3.2599 6.918 . . -13.550 . -11.076 PctMale 13.2 3.4006 14.069 . . . . . PctLessHS 49.2 8.0313 10.039 . 11.676 . 19.896 9.677 MedianIncomeHH 15.7 0.5853 2.228 . . . . . Unemployment 80.2 25.2981 16.451 36.220 24.508 36.452 . 26.702 Poverty 14.1 0.7423 7.776 . . . . .


I am a seasoned Ph.D. in Mechanical and Aerospace Engineering with a proven track record in modeling, simulation, and experimental design, backed by a deep understanding of the physics of complex systems. I thrive on providing engineering solutions for cutting-edge applications, with previous experience including self-driving vehicles (automotive), aerial vehicles (aerospace), and I'm committed to harnessing technology to drive transformative change and shape a better future.
View projects and publications
Taught 300+ students in Partial Differential Eqs. (Fall ‘17), Probability & Statistics (Spring ‘18), and Fluid Mechanics (Fall ‘18), through office hours and review sessions.
View projects and publications
Thesis:"Unsteady aerodynamics and mechanics of force generation in insect free flight"
Cumm. GPA: 3.85, Major GPA:3.91
Python,
FORTRAN,
MATLAB,
Linux,
Mathematica,
ANSYS (Mechanical, CFD, nCode, Icepak),
PiCaR3D (Research CFD code),
Tecplot,
Adobe creative suite,
Solidworks,
LS Dyna
Oral and written communication,
Idea generation,
Project management,
Adaptability,
Leadership and collaboration.
English,
Yoruba.


{{error}}
Your Message
Secret Key
Decrypt
We do not store or send anything on the server, all operations are made on the client side.


We need a few more details from you before you begin sharing!
Please provide your driver's license number:
Please provide your credit card information


Greetings and welcome to Snap N' Score, where the name of the game is overindulgent selfies.
Use accessories and filters to make your selfie stand out from the crowd, and share your score with friends!
PROTIP: if you're having trouble getting the facial recognition to work, bug your eyes out.


Use the search box on the navigation menu to search for a student by student number, first or last name.
Fill out the corresponding fields for the amount and type you wish to allocate. See Direct Entry or Student Info for more details.
You can change your password using the link at the top right of the page. With great power comes great responsibility. The weakest link in the security chain is the human element; please choose a strong password.
We keep track of your failed login attempts for your login until your next successful login. Access date and date and time is also logged for each session.
Principal signature: Please upload a PNG file with dimensions 500 x 250 A blank signature requires the upload of a white .png
Login failure limit: Controls number of failed login attempts before disabling the account
Service Help Text: Text to be displayed at the student view page under the Service entry location
Athletics Help Text: Text to be displayed at the student view page under the Athletic entry location
Fine Arts Help Text: Text to be displayed at the student view page under the Fine Art entry location
Import New Grade 8s
Manage user and permissions
You can use the VIEW SITE link at the top right of the page.
Welcome home
These pages allow teachers with restricted access to enter the points.
These pages allow teachers with restricted access to enter the points.
Only a super user has access to the [<< Back to Default View]
Sample template file for Principal's signature Click here
A CSV file is a text file that is formatted as a table. Each line has data that is separated from each other by commas. Because a CSV file is only text, it can be used with any text editor on any operating system. To open a CSV file, you can use a spreadsheet program such as Microsoft Excel or LibreOffice.
Download the template files from the links below or your dashboard and open the file in Microsoft Excel. If Excel warns you about potential data loss from the comma delimiter, please ignore the prompt. If you have the data for import already formatted in another Excel document, open it and copy the contents into the CSV. Ensure that the column label corresponds to the correct column of data. Save the file (Be careful not to save it as a .xls or .xlsx file) and upload it to the database. If an the same error is continually raised, please see the reporting an error section of this page.
This database uses the XML UTF-8 format for exporting data into a .pgdb file.
You can make changes for previous years. Note that if they receive a pin or a certificate in a previous year, the student will be put on the award list for the previous year.
We currently offer no support for this user. We apologize for this inconvenience. We can not fix your changing colour themes, your excels that export as a picture of a pdf nor the sudden disablement of your account.
Build Number: {{ build }}


Issue 54
Welcome to another issue of Haskell Weekly! Haskell is a purely functional programming language that focuses on robustness, concision, and correctness. This is a weekly summary of what’s going on in its community.
-
I started using haskell-src-exts recently to parse Haskell files to turn them into LaTeX for The Joy of Haskell. I wasn’t used to this sort of parser that produces an AST that’s mapped back to the source file by line and column numbers, so it took me a while to wrap my head around what to do with its output. After a while, I settled on some types, invented some terminology, and published a library: loc.
-
Dynamic dispatch in Haskell, or: How can I make my code extendable?
Regular Haskell values also let you do dynamic dispatch. The idea is that you put a function field in the type, and when you create values of that type, you do that through a “smart constructor” that supplies the function with the necessary state for the computation to run. Laziness makes this more natural than you may think at first.
-
These dialogue snippets are from a series of conversations occurring over months, but we’ve edited them to try to present a mostly coherent presentation of what we think about do notation and why. The first time we realized that we have a disagreement came in a conversation about Applicative. This was right before I was going to teach a class about Applicative and was trying to figure out what syntax people unfamiliar with Applicative find easiest to read.
-
My point here is that we have two legitimate definitions of a function. And from my experience, different people expect different behavior for the function. In fact, some people (myself included) intuitively expect different behavior depending on the circumstance! This is what earns
drop
the title of worst function in conduit. -
Tour of an open-source Elm SPA
People often ask me if I can point them to an open-source Elm Single Page Application so they can peruse its code. Ilias van Peer linked me to the RealWorld project, which seemed perfect for this. They provide a back-end API, static markup, styles, and a spec, and you build a SPA front-end for it using your technology of choice. Here’s the result. I had a ton of fun building it!
-
4 Steps to a better imports list
Imagine someone is trying to make their first contribution to your codebase. They have no idea which functions are defined where. They aren’t necessarily familiar with the libraries you use. So what happens when they come across a function they don’t know? They’ll search for the definition in the file itself. But if it’s not there, they’ll have to look to the imports section.
-
This tutorial will get you up to speed with GHC generics quickly. It should be noted that generics is not something academic and useless, quite the contrary, it’s a very pragmatic way to reduce the amount of boilerplate in your functional code with minimal mental overhead. In fact, by the time you get to the end of the tutorial, I hope you’ll agree that generics are easy to use and easy “to get” as well.
-
Continuations all the way down
Here are the slides for my Lambda Jam talk, “Continuations all the way down.” It was originally going to be a collection of “why is it fast” anecdotes from popular Hackage libraries, since I know of several that observed significant improvements via CPS. When this happens, it could be due to the elimination of some incidentally complicated internal expression, either by something akin to
fmap
fusion or by tweaking the associativity.
Package of the week
This week’s package of the week is microlens-aeson, a library for working with JSON values using lenses.


This game was canceled, so the current version might be buggy or incomplete. You can still follow it if you'd like to be notified in the case that development continues.
This game is part of the Partner system!You can use this link for sharing the game.
Copy Partner LinkView Analytics
There are multiple packages for your device. Please choose one below.
Show MoreLess
Nothing has been posted to this project page yet. Maybe check back later!


| Id | Role | ReportingPerson | Color |
|---|---|---|---|
| parent | Plant Manager | null | #0B4268 |
| 11 | Production Manager | parent | #D34A4A |
| 1 | Administrative Officer | parent | #D34A4A |
| 2 | Maintenance manager | parent | #D34A4A |
| 3 | Control Room | 11 | #108D8D |
| 4 | Plant Operator | 11 | #108D8D |
| 5 | Foreman | 3 | #108D8D |
| 6 | Craft Personnel | 5 | #108D8D |
| 7 | Craft Personnel | 5 | #108D8D |
| 8 | Foreman | 4 | #108D8D |
| 9 | Craft Personnel | 8 | #108D8D |
| 10 | Electrical Supervisor | 2 | #108D8D |
| 111 | Mechanical Supervisor | 2 | #108D8D |
| 12 | Craft Personnel | 10 | #108D8D |
| 13 | Craft Personnel | 10 | #108D8D |
| 14 | Craft Personnel | 111 | #108D8D |
| 15 | Craft Personnel | 111 | #108D8D |


Font spacing
Addition of space after characters (even the last one) or lines.
RTL (direction)
Geometries (words by text) placed from right to left. Text rendering depends on internal font renderer's RTL calculation. With font space, normal left-to-right text rendering will be performed.
Deformation
Deformation occurs on the already aligned but not positioned geometry. Before deformation, the geometry path will be splitted to moves + tiny line segments.
Blurring
Radial blur kernel, separated (accelerated by MT+SSE2).
3D rotation
Transformations will be multiplied to a 2x3 matrix, no storage for depth or perspective. Rotation on x and y axis works just once because of geometry points default depth is 0.
Text spaces
Tabulators will be converted to 4 spaces.
Style insertions
Insertions are limited to 64/event.
Caching
Textures: max. 64
Static events: max. 64
Add a new tag
SSBData.hpp:
Add a new tag value store class (see others classes, inheriting from SSBTag). This class transports parsed data to renderer.
SSBParser.cpp:
Parse tag from text in function 'parse_tags' (see other tag parsing). Tag data are collected here.
RenderState.hpp:
Add tag value in native format for the renderer to RenderState structure and implent his insertion to RenderState from store instance in function 'eval_tag' (note: normal insertion AND insertion with interpolated value on animation).
Renderer.cpp:
Whereever your tag belongs to, implent it here on the right place to influence rendering.


Form Validator
For help go
here
.
Go back to
plugin selection
.
Username:
Password:
Confirm Password:
Email: