Stat Student

Standard deviation (SD) is a measure of variability of the population from which a sample was drawn. The standard error (SE) is a type if SD, and takes into account variability in a sample. This article sums it up. "We can estimate how much sample means will vary from the standard deviation of this sampling distribution, which we call the standard error (SE) of the estimate of the mean."

SE = SD/√(sample size)

Simple linear regression and multiple regression are subsets of general linear models. A typical GLM has K+1 parameters and looks something like this

y=B0+B1x1+B2x2+B3x3...BkXk+epsilon

and "linear" here means linear in the parameters (x's can be of higher order and the model can be called linear). B0 is still the slope (value when all Xi=0), but now B1---Bk are partial slopes, meaning each Bi they represent the change in y per unit change in Xi, holding all other Xj constant.

GLMs an be additive (y=B0+B1x1i+B2x2i+error) which indicates the belief that the effect of the explanatory variables are additive or GLMS can include an interaction term (y=B0+B1x1+B2x2+B3X1X2), meaning the explanatory variables affect each other in an increasing or decreasing way.

Dummy variables are sometimes used when using categorical variables, like groups. The use of dummy variables get around the following problem: if we have 4 regions and code them 1,2,3,4 then we are forcing the difference between groups 2 and 3 to be the same as the difference between 1 and 4.

The idea behind the previous post can be extended for finding the pdf of a function of 2 random variables.

Now x1 and x2 are continuous random variables with joint pdf f(x1, x2). Again A is the set where this joint pdf is >0. if we have a function y1=phi1(x1, x2) and y2=phi2(x1,x2),  then if the following three conditions are met:

phi1 and phi2 are 1-1 on A

dyi/dxj are continuous on A (i,j=1,2)

the Jacobian does not = 0

(note: "phi" is used in place of the greek symbol)

then we can find f(y1,y2) using this formula:

where the  J stands for the Jacobian. The Jacobian is the name given to the determinant of a 2x2 matrix consisting of partial derivatives.

SO it's he same idea as the 1 dimensional case, but multiplying by the absolute value of the Jacobian instead of dx/dy. The hard/tedious part here is finding the corresponding support of phi(A), which is used to find marginal pdfs. 

In addition to the distribution function method talked about before, you can use a different method.

It starts out like this: X is some random variable with a density function (pdf) f(x). Then we say the set where f(x)>0 is called "A" (this is also known as the support of x. Now if we have a function of x, say Y=phi(x), then, if the following three conditions are met,

phi(x) is 1-1 on A

dy/dx is continuous on A

dx/dy does not =0

we can find f(y) (i.e. as the title says, we are trying to find the pdf of a function of random variables, and Y is a function of the random variable X) using this formula:

The first step in finding f(y) we solve the function of X for X. For example

if X is a random variable with pdf f(x)=e^(-x) for X>0, and Y=-3x+5, solving this function for X would give us x=(y-5)/-3

Note also that when x>0,  y<5

Then we use this find dx/dy, which is (-1/3)

Finally, we can develop the pdf of y by plugging in (y-5)/-3 in for x in the pdf at the beginning (e^-x), and multiply by absolute value of (dx/dy). The support of y is y<0, as we found out in step 1

Confidence Intervals and Prediction Intervals with Least Squares Regression

One of the goals of regression is to estimate a value of y for a given value of x. In addition to a point estimate of y we can develop a range of y values, basically a confidence interval for an expected value of y.

E(y) is the parameter we are using the confidence interval to estimate, with the predicted value of y (Yi-hat) serving as the point estimate for E(y). The standard error of the CI is the sqrt(MSE)*sqrt((1/n)+((Xi-x-bar)^2)/Sxx))). We can develop a range of CIs for a range of x values, which then can form confidence bands around the regression line on the scatterplot. The typical pattern of a confidence band (shaded portion shown below) is a narrowing around the mean of x and then a fanning out on both ends. This fact taken together with the std error of the CI (specifically the (Xi-Xbar)^2 term) seems to say something like we are reasonably confident about a range of y values for an x value around the mean of x, and less confident elsewhere.

Seeing the confidence band spread out as x increases gives some intuition of why extrapolation can be dangerous with linear regression, as a really wide confidence interval is basically no good.

A prediction interval (band) is a similar idea to a confidence interval (band), and is used to estimate intervals on variables. PI are used to tell you where you can expect to see the next data point sampled. A PI says if you develop many prediction intervals, you'd expect that next value to lie within that prediction interval in 95% of the samples. Since we are predicting a value of a random variable, a Prediction band is always wider than a Confidence band (dotted line in scatterplot above). The formula is basically the same, but with a 1 added under the sqrt in the std error term.

Increasing sample size, alpha level (confidence coefficient) or the range of x-values sampled will all decrease both the CI and PI.

Examining Lack of Fit with Least Squares Regression

Obviously when you develop a regression line, you want it to be a good "fit" to the data. You want the regression to explain the relationship between the independent and dependent variable (without over-fitting).

One way to examine lack of fit is to check out a residual plot, which is done after developing a regression model. A residual is simply the observed value of y minus the predicted value of y at given value of the x-variable. A scatterplot of residuals vs y-hat shows us a picture we can use to evaluate the fit of a model. Since one of the assumptions we made when first developing a regression model is that we assume the errors (residuals) have equal variance, examining a residual plot allows us to check this. 

A residual plot that would indicate the assumption of homogeneous variance has not been violated looks something like this:

There is no definite pattern here. If there were a pattern, like the residuals spread out as the values of Yi-hat increased or if the residuals looked like a parabola, the model is not a good fit, and indicates some other model (weighted or a model of higher order) would need to be used instead.

A formal test for lack of fit involves partitioning SSE (variance not explained by the model) into 2 pieces: pure error variability (SSPE) and sum of squares due to lack of fit (SSLOF). SSPE are the errors due to natural variability. In order to do this test you need at least one x value to have more than 1 observation. Writing the formula out here would be too tedious, but the basic idea here is that for each value of x where you have more than one observation:

take the mean of the associated y values

subtract the mean from the observation, square it, and sum them all, this is SSPE

 SSE-SSPE=SSLOF

divide SSPE and MSLOF by their respective degrees of freedom to get mean square pure error (MSPE) and mean square lack of fit (MSLOF)

df for SSPE is the sum of the number of y values for which you have more than 1 x observation minus 1. So if you have 5 x-values each with 2 observations, the df would be (2-1)+(2-1)+(2-1)+(2-1)+(2-1)=5

df for SSLOF is (n-2)-df(SSPE)

In the previous post, I mentioned some properties of the correlation coefficient, r. Namely, when r>0 there is a positive linear relationship and when r<0, there is a negative linear relationship, and when r=1 there is a perfect linear relationship.

Looking at the slope of the regression line can tell you basically the same thing. A positive slope tells you there is a positive relationship between X and Y and so on. So why use r?

We need this coefficient because its possible to have variables with the same slope, but have drastically different looking scatterplots, and therefore very different variation.

For example, these 2 scatterplots could have roughly the same slope but the 2nd one has much less variation, as shown by the "tightness" of the points (imagine a OLS regression line in both, the points would be much close to the line in the 2nd one).

In addition to assessing SSR and SSE there are other methods of assessing model fit: MSE and R^2.

Mean Square Error (MSE) is used to estimator of the true variance, Sigma^2. 

MSE=SSE/(n-2).

The square root of MSE is called the standard error of the estimate and is used more in the inference procedures discussed below.

R^2 (coefficient of determination) is a measure of how much of the total variation can be explained by the model.

R^2=SSR/SST

It is the percentage of the variation in Y that can be explained by the regression of Y on X. In general, an R^2 value over .8 is considered a good regression fit, while anything between .6 and .8 is only considered a moderately good fit.

Related to R^2 is the Pearson correlation coefficient, which is simply "r". r measures the strength of the linear relationship between X and Y. Linear is stressed because 2 variables could have r=0 but still be related in some non-linear way

r=B1hat*sqrt(Sxx/Sxy).

r is always between -1 and 1, and a positive r implies a positive linear relationship, while r<0 implies a negative linear relationship. r=0 implies no linear relationship and r=+/-1 implies a perfect + or - linear relationship.  r is used to develop a confidence interval for the parameter, rho, which is a measure of the true underlying correlation. We also use r to test hypotheses about the true correlation.

We also want to test the significance of our estimate of slope (rarely the intercept) and develop a confidence interval for the true value of the slope.

A 100(1-a)% CI for B1=B1hat +/- t(a/2)*(std. error of B1), where this standard error = sqrt(MSE)/sqrt(Sxx).

Another interesting inference procedure is to test if our estimate of slope, B1hat, is different from 0. If it isn't, then we might as well drop it from our study and pick a better explanatory variable because its contributing nothing. This is more important for multiple regression. 

My last post talked about variability in the least squares regression model. Specifically how total variability (SST) is composed of the variability not explained by the model (SSE) plus the variability explained by the model (SSR). But how does the model explain variability? Yes there is the vertical distance between the predicted values and mean values, but how does this distance really explain anything?

Regression analysis is a type of statistical modeling where we try and develop an equation/line based on some data we have that describes the relationship between a dependent and independent variable (more than 1 explanatory/independent variable for multiple regression). Like most things in statistics, we don't know the true underlying relationship between the data, so we estimate a line based on a sample, and account for some error in our model.

Aside from the assumption that the relationship between the variables is linear, the other basic assumptions are:

the expected value of the errors is 0

the errors have the same variance

the errors are independent of each other

the errors are normally distributed

Least squares regression is a popular type of regression (though not the only kind there is). The goal of the regression model here is to minimize the sum of squared residuals. A residual is the difference between the actual value and the predicted value: ei=(Yi-Yhat). This can be represented by a vertical distance between each observed value and the prediction line. This type of regression picks the best line, and best in this case means the line that has the shortest distance for every point. This distance is measured and squared for each point. (this distance is squared so the summation of residuals makes sense when dealing with negative values). I guess the thinking behind this is that since no model is perfect you will always have some residuals in your model, why not try to make the model as close as you can. More on this further down.

The predicted or "fitted" least squares regression model is Yhat=B0hat+B1hatX ("hat" means estimate). Again this estimated model estimates the theoretical regression line Y=B0+B1X+error. Since the expected value of the errors is assumed to be 0, we don't worry about them.

B1=Sxy/Sxx, where Sxy Sum[(Xi-Xbar)(Yi-Ybar)] and Sxx=Sum[(Xi-Xbar)^2]

Once you know B1, B0=Ybar-B1*Xbar

B0 is the intercept. This may or may not make sense to interpret, depends on the nature of the problem. B1 is the slope, it is key in allowing you to predict a value of Y given some value of X.

Things that can mess with your regression model: high leverage points and high influence points. A high leverage point is a data point that seems to be an outlier in the x direction. you can check this using a scatterplot of your data. This has the potential to alter the slope of the regression line. A high influence point does affect the slope, and is usually an outlier in both the X and Y direction.

All this concern about the slope os because slope is king in the land of regression. The slope of your regression model allows you to predict a value for Y (dependent variable) for some value of x (independent variable), which is one of the whole points of regression.

So now we have a regression model, but is it any good? we need to asses its fit. We do this by examining different types of variability:

SST - total variability

SSE - Variability due to error, not explained by the model

SSR - Variability explained by the model

SST=SSE+SSR. In words this means the total variability=variability not explained by the model+variability explained by the model. The goal is to have your model explain as much of the variability as possible (be a high % of SST) by choosing a set of independent variables while still remaining as simple as possible. It is possible to fit the model really well (have the sum of squared residuals really low), and explaining a lot of variability. This can be down by transforming the data through various methods. But the downside of that is your model will perform really badly with another sample. The model is overfitted.

Once of the things stressed in class is that regression is good for explaining the data you have at hand, but is not so good at predicting future values. This is because the nature of the relationship between variables is assumed to be static by the regression model. The relationship may change in the future ( or may not even be captured correctly in the first place) so using the model to extrapolate can be misleading.

We started class by going over our professor's Theory I final, which was a pretty good review and reminder of where I need to brush up (properties of variance, properties of covariance). Also learned his section did Ch.6 which covers transformations of random variables. Not sure if he's going to re-teach this or expect us to catch up on our own. No class next week because of MLK day, so this should give me time to read the chapter and hopefully start to understand the material before the next class.

Ch.6 in our book goes over finding the probability distribution of a function of random variables (aka transformation of random variables).The idea behind wanting to know the distribution of a function of random variables is as follows:

Ybar is an estimate of mu, but how good of an estimate?

The quantities used to estimate population parameters or to make decisions about a population are functions of the n random observations that appear in a sample

Ybar is a random variable and we cannot be certain that the error of estimation will be less than a specific value, say, B.

But if we can determine the probability distribution of the estimator Y this probability distribution can be used to determine the probability that the error of estimation is less than or equal to B.

To determine the probability distribution for a function of n random variables, Y1 , Y2 , . . . , Yn , we must find the joint probability distribution for the random variables themselves

In class we briefly went over the Distribution Function Method, which summarized is:

If Y has probability density function f(y) and if U is some function of Y, then we can find FU (u) = P(U ≤ u) directly by integrating f (y) over the region for which U ≤ u. The probability density function for U is found by differentiating FU (u). 

reading through the examples in the book the process seems to be that usually youre given some pdf like f(y) = 2y, 0<y<1 and then some other function of y like U=3y-1 and asked to find the pdf of U. So first you'd start by stating the definition: 

FU(u)=P(U≤u)

then substitute for U to get:

FU(u)=P(U≤u)=P(3Y-1≤u)

then solve for Y to get:

FU(u)=P(U≤u)=P(3Y-1≤u)=P(Y≤ ((u+1)/3))

I saw a post on this blog about keeping a statistics diary, where you observe things in the news and see how statistics can apply. I really like this idea, it sort of fits under the scope of why I wanted to start this blog.

The tension between the NYPD and mayor of NYC is a huge and complex story which I won't get into here, but what interests me is this recent NY Times OpEd which notes the NYPD are writing less tickets, making fewer arrests and basically shunning their responsibilities as police officers. It makes me wonder- are we in the middle of some weird natural experiment in which the presence and effectiveness of the police in one of the busiest and biggest cities in the world is drastically reduced? What effect, if any, will this have on the crime rate? How could we measure this? How long would any change in the crime rate take to occur?

The "if any" part is what really interests me. It would be very interesting if we observed no real difference in the crime rate, despite "less" police. It has been 2 weeks now, and there hasn't been any major news story about the skyrocketing crime rate. NYC hasn't devolved into chaos. This sort of thing seems like it could be really important for cities all over the place. Could the crime rate be kept in check with less police? If so, that's a lot of tax payer money that could be avoided or reallocated elsewhere.

It's implied (I think) that the amount of crime in a city drives the number officers needed, the level of police presence (and the police department budget). How the crime rate reacts to less policing could give some insight to this. For example if the crime says the same, it may indicate some solid evidence that the relationship between crime rate and police may need to be re-examined.

At the end of the semester I was finishing up The Signal and the Noise and was sort of caught off guard with how quickly hypothesis testing was dismissed in favor of Bayesian methods. Hypothesis testing was something we spent 2/3 of the semester on in Methods I, and here it was being reduced to eh, it's ok. So between semesters I've been trying to do a little extra reading on Bayesian methods and am finding there are broadly 2 schools of thought- frequentists and Bayesians. I sort of can't believe we weren't really introduced to the critiques of hypothesis testing in class.

This post Why continue to teach and use hypothesis testing?, the title of which made me feel I just wasted a whole bunch of time and money,  is a response to a discussion I didn't even was being had. The post and some of the comments make it seem like HT is taught because "its a thing" in statistics and its always been taught. I found this set of slides interesting, especially where it says basically frequentists won out because Bayesian analysis requires a lot of computation power, which didn't exist 70 years ago, and people at the time wanted "cookbook" procedures."

The slides go on to summarize the differences between the 2 methods. In frequentists methods:

The parameters of interest are fixed and unchanging under all realistic circumstances.

No information prior to the model specification

Statistical results assume that data were from a controlled experiment.

Nothing is more important than repeatability, no matter what we pay for it.

And Bayesian methods:

Views the world probabilistically, rather than as a set of fixed phenomena that are either known or unknown.

States prior information abounds and it is important and helpful to use it.

Are very careful about stipulating assumptions and are willing to defend them.

Believes every statistical model ever created in the history of the human race is subjective; they are willing to admit it.

Most posts on the topic conclude that statisticians should really use the best method for the particular situation and not blindly follow one method of the other all the time. I get the sense that I've barely scratched the surface of Bayesian methods, I get the general idea and know Bayes Formula, but that's about it.

This was (and still is) one of the biggest questions left over in my brain from my first semester. In Theory I we were introduced to moment generating functions and told about their uses, but the way it was taught and the fact we weren't tested on it means it's still an underdeveloped and fuzzy concept in my head. After reading around online, I think this is a pretty useful description:

The function f(X)=(X-z)n is a particularly common and useful one for various values of z and n, and in order to specify what the values of z and n take, we have to use phrases like "the nth moment of X about z", or "the nth central moment of X" which is shorthand for "the nth moment of X about E[X]", or "the nth non-central moment of X" which is shorthand for "the nth moment of X about 0".

Another really good explanation of the topic is one I got from some Penn State online course material. Basically, moment is a synonym for expected value of a random variable. You can find a lot of expected values (E(Y), E(Y2), E(Y3),E(Y4)...), and some expected values are more useful than others. For example E(Y) is the mean of a random variable and E(Y2) is useful for finding variance of a random variable:

mean = E(Y) = μ

Var(Y)  = σ2 = E(Y2) − μ2 

The mean and variance are functions of moments, and sometimes they can be hard to calculate, so instead moment generating functions (mgf)  are used. 

Let Y be a discrete random variable with probability function f(y), the the mgf of Y is: m(t) = E(etY)

If a moment-generating function exists for a random variable X, then:

The mean of Y can be found by evaluating the first derivative of the moment-generating function at t = 0. That is: μ=E(X)=M′(0)

The variance of Y can be found by evaluating the first and second derivatives of the moment-generating function at t = 0. That is: σ2=E(Y2)−[E(X)]2 =M″(0)−[M′(0)]2

So to find the mean, find the mgf, etYp(y), take the first derivative and evaluate at t=0. For the variance, again find the mgf, and evaluate the 2nd derivative at t=0 and subtract the square of the first derivative evaluated at t=0.

MGFs are useful because they allow you to prove that a random variable possesses a particular probability distribution, or that 2 random variables have the same distribution. In many cases where you need to prove that two distributions are equal, it is much easier to prove equality of the moment generating functions than to prove equality of the distribution functions. From my textbook:

"If the moment-generating functions for two random variables Y and Z are equal (for all |t| < b for some b > 0), then Y and Z must have the same probability distribution. It follows that, if we can recognize the moment-generating function of a random variable Y to be one associated with a specific distribution, then Y must have that distribution."

After this little review, it now makes sense why we would care about moments and moment generating functions. Without really knowing it, we were dealing with moments every time we calculated the mean or variance (which was a lot). As for why we care if 2 random variables have the same distribution, I turned to wikipedia. 

A sequence or other collection of random variables is independent and identically distributed (i.i.d.) if each random variable has the same probability distribution as the others and all are mutually independent. The abbreviation i.i.d. is particularly common in statistics (often as iid, sometimes written IID), where observations in a sample are often assumed to be effectively i.i.d. for the purposes of statistical inference

also from the same wikipedia page this caught my eye:

One of the simplest statistical tests, the z-test, is used to test hypotheses about means of random variables. When using the z-test, one assumes (requires) that all observations are i.i.d. in order to satisfy the conditions of the central limit theorem. 

Now that I'm between semesters I have some time to myself to organize my thoughts, re-focus and think back about what I learned over the course of my first semester of grad school. Overall though I was pretty stressed out and always pressed for time I enjoyed my classes but really disliked the textbooks used. I found myself either shunning entire sections and turning to youtube or other sites for help frequently.

The intro courses covered a lot. There was a lot of memorizing and note-card making and just a wide range of stuff to take in. Off the top of my head what I remember from Stat Theory I:

Lots of different ways to look at, think about and calculate probability

Bayes theorem

Discrete probability distributions and their means and variances (Binomial, Negative Binomial, Geometric, Hypergeometric, Poisson)

Continuous probability distributions and their means and variances  (Uniform, Normal, Gamma, Chi-Square, Exponential, Beta)

Discrete and continuous joint (multivariate) probability distributions and how to calculate their means and expected values

Covariance, correlation coefficient, conditional expectation and conditional variance

Things I wished we covered in more depth from Stat Theory I:

Moment generating functions - why are they so important? are they ever used in real life?

Conditional expectation and conditional variance - these concepts felt rushed at the end of the semester, I should probably review these before the next semester starts.

Stat Methods I was a more enjoyable class, less dry than theory I and overall more interesting. Things I remember from Methods I:

Statistic vs parameter, sampling methods & biases, describing data (distributions, histograms, box plot...) measures of central tendency and spread, more probability, CLT, sampling distributions

confidence intervals for mu and variance, hypothesis testing for mu and variance, estimating mu sigma known vs sigma unknown, comparing 2 estimates of mu, comparing variances, calculating sample size and power, interpreting p-values, estimating true proportions, comparing 2 proportions, goodness of fit tests, tests for independence and homogeneity, brief intro to regression, plus a ton more that i cant remember at the moment.

Things I wished we covered in more depth from Stat Methods I:

Thinking back i'm still a little fuzzy on the topic of statistical power. It is the probability you correctly reject a null hypothesis when it is false, and is inversely related to Beta, but it seems like a somewhat deep concept and one I should probably look into more.

How hypothesis testing relates to frequentist vs bayesian statistics. In class hypothesis testing is represented as if it is the only way to "do" statistics. Though from my reading it seems to me there is another way of thinking about things, namely Bayesian, and I wish some of the criticisms of hypothesis testing/frequentist methods was discussed more.