text:
# Lifts and Projections of Sets defined by Polynomial Inequalities

My research is on using design optimisation studies using ridge functions of quantities of interest. Ridge functions allow us to create surrogate models in lower dimensional subspaces using fewer function evaluations than would be required in the full $$n$$-dimensional space (defined to lie in the $$n$$-dimensional hypercube i.e. $$\mathbf{x} \in [-1,1]^n$$).

So for each quantity of interest, we have a ridge function $$f(\mathbf{x})\approx g(\mathbf{u})$$, where $$\mathbf{u} = U^T\mathbf{x}$$ ($$U \in \mathbb{R}^{n \times d}$$ and $$d) are the active variables in which the function varies the most. We also have the inactive variables $$\mathbf{v} = V^T\mathbf{x}$$ with which the quantity of interest does not vary a great deal, where $$V=\mathrm{null}(U^T)$$. Note, concatenating these matrices gives the orthogonal matrix $$W=[U,V]$$ and given the active and inactive variables, we can return to $$n$$-dimensional space using $$\mathbf{x}=U\mathbf{u} + V\mathbf{v}$$.

Given a ridge function $$f_1 \approx g_1(\mathbf{u}_1)$$, I can define a feasible set, say $$\Omega_1$$, in terms of its active variables by a set of $$m$$ inequalities i.e. $$\Omega_1 = \{ \mathbf{u}_1 | h_i(\mathbf{u}_1)\leq b_i, \quad i = 1,\dots,m \}$$. However, these active and inactive variables are determined by their effect on their respective quantities of interest, and so reside in different subspaces.

I am interested in what I can say about this set in the subspace defined by a different ridge function $$f_2 \approx g_2(\mathbf{u}_2)$$. That is, what is there to be said about the set $$\Omega_2 = \{ \mathbf{u}_2 | \mathbf{u}_1 \in \Omega_1 \}$$?

In the case where $$g_1(\mathbf{u}_1)$$ is a linear ridge function, the set $$\Omega_1$$ is an $$n$$-dimensional polytope (including the bounds from $$\mathbf{x}$$ being contained in the unit hypercube). This means we can use a number of available polytope projection algorithms to find the set $$\Omega_2$$ (relatively) easily.

However, I am unsure of where to even begin (or indeed if such a representation exists) in the case where $$g_1(\mathbf{u}_1)$$ is nonlinear (e.g. a quadratic). If anyone has any suggestions as to pertinent research or reading that I may look into to begin investigating this question, I would really appreciate it.

To give a visual aide, I used a very naive and inefficient sampling technique to investigate for relatively low dimensions to get some images of the sets I would expect. Feasible set defined in space $$\mathbf{u}_1$$ projected onto the space $$\mathbf{u}_2$$ The grey sets being the feasible sets and the black dashed line being the boundary of the active variables such that the condition that $$\mathbf{x}$$ being contained in the unit hypercube is met.

Thank you for taking the time to read my post.


text:
# Math Help - how fast distance is changing

1. ## how fast distance is changing

A plane flies directly over an observation point. The plane has constant altitude of 5000 ft. How fast is the distance between the plane and the point changing, at the instant the plane is 600ft past the overhead position and the plane is traveling 300 ft miles per hour?

NOTE: 5000ft=0.947 mile, and 600ft=0.1136 mile

2. Use Pythagoras.

$D^{2}=x^{2}+y^{2}$

Differentiate wrt time.

$D\cdot\frac{dD}{dt}=x\cdot\frac{dx}{dt}+y\cdot\fra c{dy}{dt}$

Let D be the hypoteneuse formed by the triangle we get.

You can find D by just using Pythagoras right out given x=600 and y=5000.

Remember, consistent units. That is the trick here. Be careful.

Fortunately, you're given the conversions.

Also y is constant, so dy/dt=0

You have all the knowns. Enter them in and solve for dD/dt

3. i am still a little confused....can you help me set this up a little further? im not sure where 300 mph comes in

4. 300 is dx/dt. That is the rate of the plane flying horizontally. The rate of change of x.

format: plain


text:
# If y varies inversely as x, how do you find the constant of variation and the inverse variation equation given y=90 when x =6?

$y = \frac{540}{x} \mathmr{and} x = \frac{540}{y}$.
$y = \frac{k}{x}$, where k is the constant ox proportionality.
In this unique relation, the function ${y}^{- 1}$ is same as the function y.


text:
Fitting ordered probit models with endogenous covariates with Stata’s gsem command

The new command gsem allows us to fit a wide variety of models; among the many possibilities, we can account for endogeneity on different models. As an example, I will fit an ordinal model with endogenous covariates.

Parameterizations for an ordinal probit model

The ordinal probit model is used to model ordinal dependent variables. In the usual parameterization, we assume that there is an underlying linear regression, which relates an unobserved continuous variable $$y^*$$ to the covariates $$x$$.

$y^*_{i} = x_{i}\gamma + u_i$

The observed dependent variable $$y$$ relates to $$y^*$$ through a series of cut-points $$-\infty =\kappa_0<\kappa_1<\dots< \kappa_m=+\infty$$ , as follows:

$y_{i} = j {\mbox{ if }} \kappa_{j-1} < y^*_{i} \leq \kappa_j$

Provided that the variance of $$u_i$$ can’t be identified from the observed data, it is assumed to be equal to one. However, we can consider a re-scaled parameterization for the same model; a straightforward way of seeing this, is by noting that, for any positive number $$M$$:

$\kappa_{j-1} < y^*_{i} \leq \kappa_j \iff M\kappa_{j-1} < M y^*_{i} \leq M\kappa_j$

that is,

$\kappa_{j-1} < x_i\gamma + u_i \leq \kappa_j \iff M\kappa_{j-1}< x_i(M\gamma) + Mu_i \leq M\kappa_j$

In other words, if the model is identified, it can be represented by multiplying the unobserved variable $$y$$ by a positive number, and this will mean that the standard error of the residual component, the coefficients, and the cut-points will be multiplied by this number.

Let me show you an example; I will first fit a standard ordinal probit model, both with oprobit and with gsem. Then, I will use gsem to fit an ordinal probit model where the residual term for the underlying linear regression has a standard deviation equal to 2. I will do this by introducing a latent variable $$L$$, with variance 1, and coefficient $$\sqrt 3$$. This will be added to the underlying latent residual, with variance 1; then, the ‘new’ residual term will have variance equal to $$1+((\sqrt 3)^2\times Var(L))= 4$$, so the standard deviation will be 2. We will see that as a result, the coefficients, as well as the cut-points, will be multiplied by 2.

. sysuse auto, clear
(1978 Automobile Data)

. oprobit rep mpg disp , nolog

Ordered probit regression                         Number of obs   =         69
LR chi2(2)      =      14.68
Prob > chi2     =     0.0006
Log likelihood = -86.352646                       Pseudo R2       =     0.0783

------------------------------------------------------------------------------
rep78 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
mpg |   .0497185   .0355452     1.40   0.162    -.0199487    .1193858
displacement |  -.0029884   .0021498    -1.39   0.165     -.007202    .0012252
-------------+----------------------------------------------------------------
/cut1 |  -1.570496   1.146391                      -3.81738    .6763888
/cut2 |  -.7295982   1.122361                     -2.929386     1.47019
/cut3 |   .6580529   1.107838                     -1.513269    2.829375
/cut4 |    1.60884   1.117905                     -.5822132    3.799892
------------------------------------------------------------------------------

. gsem (rep <- mpg disp, oprobit), nolog

Generalized structural equation model             Number of obs   =         69
Log likelihood = -86.352646

--------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
rep78 <-       |
mpg |   .0497185   .0355452     1.40   0.162    -.0199487    .1193858
displacement |  -.0029884   .0021498    -1.39   0.165     -.007202    .0012252
---------------+----------------------------------------------------------------
rep78          |
/cut1 |  -1.570496   1.146391    -1.37   0.171     -3.81738    .6763888
/cut2 |  -.7295982   1.122361    -0.65   0.516    -2.929386     1.47019
/cut3 |   .6580529   1.107838     0.59   0.553    -1.513269    2.829375
/cut4 |    1.60884   1.117905     1.44   0.150    -.5822132    3.799892
--------------------------------------------------------------------------------

. local a = sqrt(3)

. gsem (rep <- mpg disp L@a'), oprobit var(L@1) nolog

Generalized structural equation model             Number of obs   =         69
Log likelihood = -86.353008

( 1)  [rep78]L = 1.732051
( 2)  [var(L)]_cons = 1
--------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
rep78 <-       |
mpg |    .099532     .07113     1.40   0.162    -.0398802    .2389442
displacement |  -.0059739   .0043002    -1.39   0.165    -.0144022    .0024544
L |   1.732051  (constrained)
---------------+----------------------------------------------------------------
rep78          |
/cut1 |  -3.138491   2.293613    -1.37   0.171     -7.63389    1.356907
/cut2 |  -1.456712   2.245565    -0.65   0.517    -5.857938    2.944513
/cut3 |   1.318568    2.21653     0.59   0.552     -3.02575    5.662887
/cut4 |   3.220004   2.236599     1.44   0.150     -1.16365    7.603657
---------------+----------------------------------------------------------------
var(L)|          1  (constrained)
--------------------------------------------------------------------------------

Ordinal probit model with endogenous covariates

This model is defined analogously to the model fitted by -ivprobit- for probit models with endogenous covariates; we assume Read more…

Categories: Statistics Tags:

Export tables to Excel

There is a new command in Stata 13, putexcel, that allows you to easily export matrices, expressions, and stored results to an Excel file. Combining putexcel with a Stata command’s stored results allows you to create the table displayed in your Stata Results window in an Excel file.

A stored result is simply a scalar, macro, or matrix stored in memory after you run a Stata command. The two main types of stored results are e-class (for estimation commands) and r-class (for general commands). You can list a command’s stored results after it has been run by typing ereturn list (for estimation commands) and return list (for general commands). Let’s try a simple example by loading the auto dataset and running correlate on the variables foreign and mpg

. sysuse auto
(1978 Automobile Data)

. correlate foreign mpg
(obs=74)

|  foreign      mpg
-------------+------------------
foreign |   1.0000
mpg |   0.3934   1.0000

Because correlate is not an estimation command, use the return list command to see its stored results.

. return list

scalars:
r(N) =  74
r(rho) =  .3933974152205484

matrices:
r(C) :  2 x 2

Now we can use putexcel to export these results to Excel. The basic syntax of putexcel is

putexcel excel_cell=(expression) … using filename [, options]

If you are working with matrices, the syntax is

putexcel excel_cell=matrix(expression) … using filename [, options]

It is easy to build the above syntax in the putexcel dialog. There is a helpful video on Youtube about the dialog here. Let’s list the matrix r(C) to see what it contains.

. matrix list r(C)

symmetric r(C)[2,2]
foreign        mpg
foreign          1
mpg  .39339742          1

To re-create the table in Excel, we need to export the matrix r(C) with the matrix row and column names. The command to type in your Stata Command window is

putexcel A1=matrix(r(C), names) using corr

Note that to export the matrix row and column names, Read more…

Categories: Programming Tags:

Measures of effect size in Stata 13

Today I want to talk about effect sizes such as Cohen’s d, Hedges’s g, Glass’s Δ, η2, and ω2. Effects sizes concern rescaling parameter estimates to make them easier to interpret, especially in terms of practical significance.

Many researchers in psychology and education advocate reporting of effect sizes, professional organizations such as the American Psychological Association (APA) and the American Educational Research Association (AERA) strongly recommend their reporting, and professional journals such as the Journal of Experimental Psychology: Applied and Educational and Psychological Measurement require that they be reported.

Anyway, today I want to show you

1. What effect sizes are.
2. How to calculate effect sizes and their confidence intervals in Stata.
3. How to calculate bootstrap confidence intervals for those effect sizes.
4. How to use Stata’s effect-size calculator.

1. What are effect sizes?

The importance of research results is often assessed by statistical significance, usually that the p-value is less than 0.05. P-values and statistical significance, however, don’t tell us anything about practical significance.

What if I told you that I had developed a new weight-loss pill and that the difference between the average weight loss for people who took the pill and the those who took a placebo was statistically significant? Would you buy my new pill? If you were overweight, you might reply, “Of course! I’ll take two bottles and a large order of french fries to go!”. Now let me add that the average difference in weight loss was only one pound over the year. Still interested? My results may be statistically significant but they are not practically significant.

Or what if I told you that the difference in weight loss was not statistically significant — the p-value was “only” 0.06 — but the average difference over the year was 20 pounds? You might very well be interested in that pill.

The size of the effect tells us about the practical significance. P-values do not assess practical significance.

All of which is to say, one should report parameter estimates along with statistical significance.

In my examples above, you knew that 1 pound over the year is small and 20 pounds is large because you are familiar with human weights.

In another context, 1 pound might be large, and in yet another, 20 pounds small.

Formal measures of effects sizes are thus usually presented in unit-free but easy-to-interpret form, such as standardized differences and proportions of variability explained.

The “d” family

Effect sizes that measure the scaled difference between means belong to the “d” family. The generic formula is

The estimators differ in terms of how sigma is calculated.

Cohen’s d, for instance, uses the pooled sample standard deviation.

Hedges’s g incorporates an adjustment which removes the bias of Cohen’s d.

Glass’s Δ was originally developed in the context of Read more…

Categories: Statistics Tags:

Stata 13 ships June 24

There’s a new release of Stata. You can order it now, it starts shipping on June 24, and you can find out about it at www.stata.com/stata13/.

Well, we sure haven’t made that sound exciting when, in fact, Stata 13 is a big — we mean really BIG — release, and we really do want to tell you about it.

Rather than summarizing, however, we’ll send you to the website, which in addition to the standard marketing materials, has technical sheets, demonstrations, and even videos of the new features.

And all 11,000 pages of the manuals are now online.

Update on the Stata YouTube Channel

What is it about round numbers that compels us to pause and reflect? We celebrate 20-year school reunions, 25-year wedding anniversaries, 50th birthdays and other similar milestones. I don’t know the answer but the Stata YouTube Channel recently passed several milestones – more than 1500 subscribers, over 50,000 video views and it was launched six months ago. We felt the need for a small celebration to mark the occasion, and I thought that I would give you a brief update.

I could tell you about re-recording the original 24 videos with a larger font to make them easier to read. I could tell you about the hardware and software that we use to record them including our experiments with various condenser and dynamic microphones. I could share quotes from some of the nice messages we’ve received. But I think it would be more fun to talk about….you!

YouTube collects data about the number of views each video receives as well as summary data about who, what, when, where, and how you are watching them. There is no need to be concerned about your privacy; there are no personal identifiers of any kind associated with these data. But the summary data are interesting, and I thought it might be fun to share some of the data with you.

Who’s watching?

Figure 1 shows the age distribution of Stata YouTube Channel viewers. If you have ever attended a Stata Conference, you will not be surprised by this graph…until you notice the age group at the bottom. I would not have guessed that 13-17 year olds are watching our videos. Perhaps they saw Stata in the movie “Moneyball” with Brad Pitt and wanted to learn more. Or maybe they were influenced by the latest fashion craze sweeping the youth of the world.

What are you watching?

We have posted more than 50 videos over a wide range of topics. Figure 2 shows the total number of views for the ten most popular videos. The more popular of the ten are about broad topics. These broader videos are mostly older and have thus had time to accumulate more views.

Even so, these videos receive more views per day currently than do the special topic videos that have been posted more recently. This supports my belief that Stata YouTube Channel viewers tend to be relatively new Stata users who want to learn about general topics, and that means more generic videos in the future. So you and your two post-docs will just have to read the manual if you want to learn how to fit asymmetric power ARCH models with outer-product gradient standard errors.

When are you watching?

We usually post new videos on Tuesday mornings which might lead you to believe that the peak viewing day would also be Tuesday. Figure 3, however, shows us that the average number of views per day (vpd) is higher on Wednesdays at 420 vpd and in fact peaks on Thursdays at 430 vpd before declining Friday through Sunday.

Figure 4 also shows us that late September may have been not the best time to launch the Stata YouTube Channel. Our early momentum in September and October slowed during the November and December holiday seasons. We were, however, pleased to see that 49 of you spent New Years Eve watching our videos. Perhaps next year we’ll prepare something more festive just for you!

Where are you watching?

What do the Czech Republic, Pakistan, Uganda, Madagascar, the United Kingdom, the Bahamas, the United States, Montenegro, and Italy have in common? Correct! They are all countries in which you are watching our videos. They are also locations depicted in one of my favorite action films but I’ll leave that to the trivia buffs. I think the most exciting information that we found in our data is that the Stata YouTube Channel is being viewed in 164 countries!

You might not be surprised to learn that roughly half of the people watching the videos live in the United States, the United Kingdom, or Canada. The results may be unexpected when we consider the “view rate” defined as the number of views per 100,000 residents. Figure 5 shows the top 20 countries ranked by view rate for countries with at least four million residents. Denmark had the highest view rate which was nearly twice the rate of Norway which had the second highest view rate. The view rate in Denmark was more than three times the rate in the US and the UK.

How are you watching?

You might think that I would have anything to report about “how” you are watching the videos, but it turns out that 5.2% of you are watching on mobile devices. Perhaps this explains the 13-17 year old demographic or the 49 people watching on New Year’s Eve. Or maybe we are helping you pass the time in the dentist office waiting room.

Final thoughts

Six months isn’t much of a milestone. We Stata folk will use any excuse to break out the cake and ice cream. Even so, the Stata YouTube Channel began as an experiment and often experiments do not work out as we would like. This experiment has exceeded our expectations and, as a result, we have started taking requests for videos on our Facebook page and we’ll be adding more videos every week. So thanks for watching and stay tuned!

Now if you will excuse me, I’m going to get some cake and ice cream.

Categories: Resources Tags:

Multilevel linear models in Stata, part 2: Longitudinal data

In my last posting, I introduced you to the concepts of hierarchical or “multilevel” data. In today’s post, I’d like to show you how to use multilevel modeling techniques to analyse longitudinal data with Stata’s xtmixed command.

Last time, we noticed that our data had two features. First, we noticed that the means within each level of the hierarchy were different from each other and we incorporated that into our data analysis by fitting a “variance component” model using Stata’s xtmixed command.

The second feature that we noticed is that repeated measurement of GSP showed an upward trend. We’ll pick up where we left off last time and stick to the concepts again and you can refer to the references at the end to learn more about the details.

The videos

Stata has a very friendly dialog box that can assist you in building multilevel models. If you would like a brief introduction using the GUI, you can watch a demonstration on Stata’s YouTube Channel:

Introduction to multilevel linear models in Stata, part 2: Longitudinal data

Longitudinal data

I’m often asked by beginning data analysts – “What’s the difference between longitudinal data and time-series data? Aren’t they the same thing?”.

The confusion is understandable — both types of data involve some measurement of time. But the answer is no, they are not the same thing.

Univariate time series data typically arise from the collection of many data points over time from a single source, such as from a person, country, financial instrument, etc.

Longitudinal data typically arise from collecting a few observations over time from many sources, such as a few blood pressure measurements from many people.

There are some multivariate time series that blur this distinction but a rule of thumb for distinguishing between the two is that time series have more repeated observations than subjects while longitudinal data have more subjects than repeated observations.

Because our GSP data from last time involve 17 measurements from 48 states (more sources than measurements), we will treat them as longitudinal data.

Random intercept models

As I mentioned last time, repeated observations on a group of individuals can be conceptualized as multilevel data and modeled just as any other multilevel data. We left off last time with a variance component model for GSP (Gross State Product, logged) and noted that our model assumed a constant GSP over time while the data showed a clear upward trend.

If we consider a single observation and think about our model, nothing in the fixed or random part of the models is a function of time.

Let’s begin by adding the variable year to the fixed part of our model.

As we expected, our grand mean has become a linear regression which more accurately reflects the change over time in GSP. What might be unexpected is that each state’s and region’s mean has changed as well and now has the same slope as the regression line. This is because none of the random components of our model are a function of time. Let’s fit this model with the xtmixed command:

. xtmixed gsp year, || region: || state:

------------------------------------------------------------------------------
gsp |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
year |   .0274903   .0005247    52.39   0.000     .0264618    .0285188
_cons |  -43.71617   1.067718   -40.94   0.000    -45.80886   -41.62348
------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
region: Identity             |
sd(_cons) |   .6615238   .2038949      .3615664    1.210327
-----------------------------+------------------------------------------------
state: Identity              |
sd(_cons) |   .7805107   .0885788      .6248525    .9749452
-----------------------------+------------------------------------------------
sd(Residual) |   .0734343   .0018737      .0698522    .0772001
------------------------------------------------------------------------------

The fixed part of our model now displays an estimate of the intercept (_cons = -43.7) and the slope (year = 0.027). Let’s graph the model for Region 7 and see if it fits the data better than the variance component model.

predict GrandMean, xb
label var GrandMean "GrandMean"
predict RegionEffect, reffects level(region)
predict StateEffect, reffects level(state)
gen RegionMean = GrandMean + RegionEffect
gen StateMean = GrandMean + RegionEffect + StateEffect

twoway  (line GrandMean year, lcolor(black) lwidth(thick))      ///
(line RegionMean year, lcolor(blue) lwidth(medthick))   ///
(line StateMean year, lcolor(green) connect(ascending)) ///
(scatter gsp year, mcolor(red) msize(medsmall))         ///
if region ==7,                                          ///
ytitle(log(Gross State Product), margin(medsmall))      ///
legend(cols(4) size(small))                             ///
title("Multilevel Model of GSP for Region 7", size(medsmall))

That looks like a much better fit than our variance-components model from last time. Perhaps I should leave well enough alone, but I can’t help noticing that the slopes of the green lines for each state don’t fit as well as they could. The top green line fits nicely but the second from the top looks like it slopes upward more than is necessary. That’s the best fit we can achieve if the regression lines are forced to be parallel to each other. But what if the lines were not forced to be parallel? What if we could fit a “mini-regression model” for each state within the context of my overall multilevel model. Well, good news — we can!

Random slope models

By introducing the variable year to the fixed part of the model, we turned our grand mean into a regression line. Next I’d like to incorporate the variable year into the random part of the model. By introducing a fourth random component that is a function of time, I am effectively estimating a separate regression line within each state.

Notice that the size of the new, brown deviation u1ij. is a function of time. If the observation were one year to the left, u1ij. would be smaller and if the observation were one year to the right, u1ij.would be larger.

It is common to “center” the time variable before fitting these kinds of models. Explaining why is for another day. The quick answer is that, at some point during the fitting of the model, Stata will have to compute the equivalent of the inverse of the square of year. For the year 1986 this turns out to be 2.535e-07. That’s a fairly small number and if we multiply it by another small number…well, you get the idea. By centering age (e.g. cyear = year – 1978), we get a more reasonable number for 1986 (0.01). (Hint: If you have problems with your model converging and you have large values for time, try centering them. It won’t always help, but it might).

So let’s center our year variable by subtracting 1978 and fit a model that includes a random slope.

gen cyear = year - 1978
xtmixed gsp cyear, || region: || state: cyear, cov(indep)

I’ve color-coded the output so that we can match each part of the output back to the model and the graph. The fixed part of the model appears in the top table and it looks like any other simple linear regression model. The random part of the model is definitely more complicated. If you get lost, look back at the graphic of the deviations and remind yourself that we have simply partitioned the deviation of each observation into four components. If we did this for every observation, the standard deviations in our output are simply the average of those deviations.

Let’s look at a graph of our new “random slope” model for Region 7 and see how well it fits our data.

predict GrandMean, xb
label var GrandMean "GrandMean"
predict RegionEffect, reffects level(region)
predict StateEffect_year StateEffect_cons, reffects level(state)

gen RegionMean = GrandMean + RegionEffect
gen StateMean_cons = GrandMean + RegionEffect + StateEffect_cons
gen StateMean_year = GrandMean + RegionEffect + StateEffect_cons + ///
(cyear*StateEffect_year)

twoway  (line GrandMean cyear, lcolor(black) lwidth(thick))             ///
(line RegionMean cyear, lcolor(blue) lwidth(medthick))          ///
(line StateMean_cons cyear, lcolor(green) connect(ascending))   ///
(line StateMean_year cyear, lcolor(brown) connect(ascending))   ///
(scatter gsp cyear, mcolor(red) msize(medsmall))                ///
if region ==7,                                                  ///
ytitle(log(Gross State Product), margin(medsmall))              ///
legend(cols(3) size(small))                                     ///
title("Multilevel Model of GSP for Region 7", size(medsmall))

The top brown line fits the data slightly better, but the brown line below it (second from the top) is a much better fit. Mission accomplished!

Where do we go from here?

I hope I have been able to convince you that multilevel modeling is easy using Stata’s xtmixed command and that this is a tool that you will want to add to your kit. I would love to say something like “And that’s all there is to it. Go forth and build models!”, but I would be remiss if I didn’t point out that I have glossed over many critical topics.

In our GSP example, we would still like to consider the impact of other independent variables. I haven’t mentioned choice of estimation methods (ML or REML in the case of xtmixed). I’ve assessed the fit of our models by looking at graphs, an approach important but incomplete. We haven’t thought about hypothesis testing. Oh — and, all the usual residual diagnostics for linear regression such as checking for outliers, influential observations, heteroskedasticity and normality still apply….times four! But now that you understand the concepts and some of the mechanics, it shouldn’t be difficult to fill in the details. If you’d like to learn more, check out the links below.

I hope this was helpful…thanks for stopping by.

If you’d like to learn more about modeling multilevel and longitudinal data, check out

Multilevel and Longitudinal Modeling Using Stata, Third Edition
Volume I: Continuous Responses
Volume II: Categorical Responses, Counts, and Survival
by Sophia Rabe-Hesketh and Anders Skrondal

or sign up for our popular public training course Multilevel/Mixed Models Using Stata.

Categories: Statistics Tags:

Multilevel linear models in Stata, part 1: Components of variance

In the last 15-20 years multilevel modeling has evolved from a specialty area of statistical research into a standard analytical tool used by many applied researchers.

Stata has a lot of multilevel modeling capababilities.

I want to show you how easy it is to fit multilevel models in Stata. Along the way, we’ll unavoidably introduce some of the jargon of multilevel modeling.

I’m going to focus on concepts and ignore many of the details that would be part of a formal data analysis. I’ll give you some suggestions for learning more at the end of the post.

The videos

Stata has a friendly dialog box that can assist you in building multilevel models. If you would like a brief introduction using the GUI, you can watch a demonstration on Stata’s YouTube Channel:

Introduction to multilevel linear models in Stata, part 1: The xtmixed command

Multilevel data

Multilevel data are characterized by a hierarchical structure. A classic example is children nested within classrooms and classrooms nested within schools. The test scores of students within the same classroom may be correlated due to exposure to the same teacher or textbook. Likewise, the average test scores of classes might be correlated within a school due to the similar socioeconomic level of the students.

You may have run across datasets with these kinds of structures in your own work. For our example, I would like to use a dataset that has both longitudinal and classical hierarchical features. You can access this dataset from within Stata by typing the following command:

use http://www.stata-press.com/data/r12/productivity.dta

We are going to build a model of gross state product for 48 states in the USA measured annually from 1970 to 1986. The states have been grouped into nine regions based on their economic similarity. For distributional reasons, we will be modeling the logarithm of annual Gross State Product (GSP) but in the interest of readability, I will simply refer to the dependent variable as GSP.

. describe gsp year state region

storage  display     value
variable name   type   format      label      variable label
-----------------------------------------------------------------------------
gsp             float  %9.0g                  log(gross state product)
year            int    %9.0g                  years 1970-1986
state           byte   %9.0g                  states 1-48
region          byte   %9.0g                  regions 1-9

Let’s look at a graph of these data to see what we’re working with.

twoway (line gsp year, connect(ascending)), ///
by(region, title("log(Gross State Product) by Region", size(medsmall)))

Each line represents the trajectory of a state’s (log) GSP over the years 1970 to 1986. The first thing I notice is that the groups of lines are different in each of the nine regions. Some groups of lines seem higher and some groups seem lower. The second thing that I notice is that the slopes of the lines are not the same. I’d like to incorporate those attributes of the data into my model.

Components of variance

Let’s tackle the vertical differences in the groups of lines first. If we think about the hierarchical structure of these data, I have repeated observations nested within states which are in turn nested within regions. I used color to keep track of the data hierarchy.

We could compute the mean GSP within each state and note that the observations within in each state vary about their state mean.

Likewise, we could compute the mean GSP within each region and note that the state means vary about their regional mean.

We could also compute a grand mean and note that the regional means vary about the grand mean.

Next, let’s introduce some notation to help us keep track of our mutlilevel structure. In the jargon of multilevel modelling, the repeated measurements of GSP are described as “level 1”, the states are referred to as “level 2” and the regions are “level 3”. I can add a three-part subscript to each observation to keep track of its place in the hierarchy.

Now let’s think about our model. The simplest regression model is the intercept-only model which is equivalent to the sample mean. The sample mean is the “fixed” part of the model and the difference between the observation and the mean is the residual or “random” part of the model. Econometricians often prefer the term “disturbance”. I’m going to use the symbol μ to denote the fixed part of the model. μ could represent something as simple as the sample mean or it could represent a collection of independent variables and their parameters.

Each observation can then be described in terms of its deviation from the fixed part of the model.

If we computed this deviation of each observation, we could estimate the variability of those deviations. Let’s try that for our data using Stata’s xtmixed command to fit the model:

. xtmixed gsp

Mixed-effects ML regression                     Number of obs      =       816

Wald chi2(0)       =         .
Log likelihood = -1174.4175                     Prob > chi2        =         .

------------------------------------------------------------------------------
gsp |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons |   10.50885   .0357249   294.16   0.000     10.43883    10.57887
------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
sd(Residual) |   1.020506   .0252613      .9721766    1.071238
------------------------------------------------------------------------------

The top table in the output shows the fixed part of the model which looks like any other regression output from Stata, and the bottom table displays the random part of the model. Let’s look at a graph of our model along with the raw data and interpret our results.

predict GrandMean, xb
label var GrandMean "GrandMean"
twoway  (line GrandMean year, lcolor(black) lwidth(thick))              ///
(scatter gsp year, mcolor(red) msize(tiny)),                    ///
ytitle(log(Gross State Product), margin(medsmall))              ///
legend(cols(4) size(small))                                     ///
title("GSP for 1970-1986 by Region", size(medsmall))

The thick black line in the center of the graph is the estimate of _cons, which is an estimate of the fixed part of model for GSP. In this simple model, _cons is the sample mean which is equal to 10.51. In “Random-effects Parameters” section of the output, sd(Residual) is the average vertical distance between each observation (the red dots) and fixed part of the model (the black line). In this model, sd(Residual) is the estimate of the sample standard deviation which equals 1.02.

At this point you may be thinking to yourself – “That’s not very interesting – I could have done that with Stata’s summarize command”. And you would be correct.

. summ gsp

Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
gsp |       816    10.50885    1.021132    8.37885   13.04882

But here’s where it does become interesting. Let’s make the random part of the model more complex to account for the hierarchical structure of the data. Consider a single observation, yijk and take another look at its residual.

The observation deviates from its state mean by an amount that we will denote eijk. The observation’s state mean deviates from the the regionals mean uij. and the observation’s regional mean deviates from the fixed part of the model, μ, by an amount that we will denote ui... We have partitioned the observation’s residual into three parts, aka “components”, that describe its magnitude relative to the state, region and grand means. If we calculated this set of residuals for each observation, wecould estimate the variability of those residuals and make distributional assumptions about them.

These kinds of models are often called “variance component” models because they estimate the variability accounted for by each level of the hierarchy. We can estimate a variance component model for GSP using Stata’s xtmixed command:

xtmixed gsp, || region: || state:

------------------------------------------------------------------------------
gsp |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons |   10.65961   .2503806    42.57   0.000     10.16887    11.15035
------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
region: Identity             |
sd(_cons) |   .6615227   .2038944       .361566    1.210325
-----------------------------+------------------------------------------------
state: Identity              |
sd(_cons) |   .7797837   .0886614      .6240114    .9744415
-----------------------------+------------------------------------------------
sd(Residual) |   .1570457   .0040071       .149385    .1650992
------------------------------------------------------------------------------

The fixed part of the model, _cons, is still the sample mean. But now there are three parameters estimates in the bottom table labeled “Random-effects Parameters”. Each quantifies the average deviation at each level of the hierarchy.

Let’s graph the predictions from our model and see how well they fit the data.

predict GrandMean, xb
label var GrandMean "GrandMean"
predict RegionEffect, reffects level(region)
predict StateEffect, reffects level(state)
gen RegionMean = GrandMean + RegionEffect
gen StateMean = GrandMean + RegionEffect + StateEffect

twoway  (line GrandMean year, lcolor(black) lwidth(thick))      ///
(line RegionMean year, lcolor(blue) lwidth(medthick))   ///
(line StateMean year, lcolor(green) connect(ascending)) ///
(scatter gsp year, mcolor(red) msize(tiny)),            ///
ytitle(log(Gross State Product), margin(medsmall))      ///
legend(cols(4) size(small))                             ///
by(region, title("Multilevel Model of GSP by Region", size(medsmall)))

Wow – that’s a nice graph if I do say so myself. It would be impressive for a report or publication, but it’s a little tough to read with all nine regions displayed at once. Let’s take a closer look at Region 7 instead.

twoway  (line GrandMean year, lcolor(black) lwidth(thick))      ///
(line RegionMean year, lcolor(blue) lwidth(medthick))   ///
(line StateMean year, lcolor(green) connect(ascending)) ///
(scatter gsp year, mcolor(red) msize(medsmall))         ///
if region ==7,                                          ///
ytitle(log(Gross State Product), margin(medsmall))      ///
legend(cols(4) size(small))                             ///
title("Multilevel Model of GSP for Region 7", size(medsmall))

The red dots are the observations of GSP for each state within Region 7. The green lines are the estimated mean GSP within each State and the blue line is the estimated mean GSP within Region 7. The thick black line in the center is the overall grand mean for all nine regions. The model appears to fit the data fairly well but I can’t help noticing that the red dots seem to have an upward slant to them. Our model predicts that GSP is constant within each state and region from 1970 to 1986 when clearly the data show an upward trend.

So we’ve tackled the first feature of our data. We’ve succesfully incorporated the basic hierarchical structure into our model by fitting a variance componentis using Stata’s xtmixed command. But our graph tells us that we aren’t finished yet.

Next time we’ll tackle the second feature of our data — the longitudinal nature of the observations.

If you’d like to learn more about modelling multilevel and longitudinal data, check out

Multilevel and Longitudinal Modeling Using Stata, Third Edition
Volume I: Continuous Responses
Volume II: Categorical Responses, Counts, and Survival
by Sophia Rabe-Hesketh and Anders Skrondal

or sign up for our popular public training course “Multilevel/Mixed Models Using Stata“.

There’s a course coming up in Washington, DC on February 7-8, 2013.

Categories: Statistics Tags:

Using Stata’s random-number generators, part 4, details

For those interested in how pseudo random number generators work, I just wrote something on Statalist which you can see in the Statalist archives by clicking the link even if you do not subscribe:

http://www.stata.com/statalist/archive/2012-10/msg01129.html

To remind you, I’ve been writing about how to use random-number generators in parts 1, 2, and 3, and I still have one more posting I want to write on the subject. What I just wrote on Statalist, however, is about how random-number generators work, and I think you will find it interesting.

To find out more about Statalist, see

Statalist

How to successfully ask a question on Statalist

Categories: Numerical Analysis Tags:

Using Stata’s SEM features to model the Beck Depression Inventory

I just got back from the 2012 Stata Conference in San Diego where I gave a talk on Psychometric Analysis Using Stata and from the 2012 American Psychological Association Meeting in Orlando. Stata’s structural equation modeling (SEM) builder was popular at both meetings and I wanted to show you how easy it is to use. If you are not familiar with the basics of SEM, please refer to the references at the end of the post. My goal is simply to show you how to use the SEM builder assuming that you already know something about SEM. If you would like to view a video demonstration of the SEM builder, please click the play button below:

The data used here and for the silly examples in my talk were simulated to resemble one of the most commonly used measures of depression: the Beck Depression Inventory (BDI). If you find these data too silly or not relevant to your own research, you could instead imagine it being a set of questions to measure mathematical ability, the ability to use a statistical package, or whatever you wanted.

The Beck Depression Inventory

Originally published by Aaron Beck and colleagues in 1961, the BDI marked an important change in the conceptualization of depression from a psychoanalytic perspective to a cognitive/behavioral perspective. It was also a landmark in the measurement of depression shifting from lengthy, expensive interviews with a psychiatrist to a brief, inexpensive questionnaire that could be scored and quantified. The original inventory consisted of 21 questions each allowing ordinal responses of increasing symptom severity from 0-3. The sum of the responses could then be used to classify a respondent’s depressive symptoms as none, mild, moderate or severe. Many studies have demonstrated that the BDI has good psychometric properties such as high test-retest reliability and the scores correlate well with the assessments of psychiatrists and psychologists. The 21 questions can also be grouped into two subscales. The affective scale includes questions like “I feel sad” and “I feel like a failure” that quantify emotional symptoms of depression. The somatic or physical scale includes questions like “I have lost my appetite” and “I have trouble sleeping” that quantify physical symptoms of depression. Since its original publication, the BDI has undergone two revisions in response to the American Psychiatric Association’s (APA) Diagnostic and Statistical Manuals (DSM) and the BDI-II remains very popular.

The Stata Depression Inventory

Since the BDI is a copyrighted psychometric instrument, I created a fictitious instrument called the “Stata Depression Inventory”. It consists of 20 questions each beginning with the phrase “My statistical software makes me…”. The individual questions are listed in the variable labels below.

. describe qu1-qu20

variable  storage  display    value
name       type   format     label      variable label
------------------------------------------------------------------------------
qu1         byte   %16.0g     response   ...feel sad
qu2         byte   %16.0g     response   ...feel pessimistic about the future
qu3         byte   %16.0g     response   ...feel like a failure
qu4         byte   %16.0g     response   ...feel dissatisfied
qu5         byte   %16.0g     response   ...feel guilty or unworthy
qu6         byte   %16.0g     response   ...feel that I am being punished
qu7         byte   %16.0g     response   ...feel disappointed in myself
qu8         byte   %16.0g     response   ...feel am very critical of myself
qu9         byte   %16.0g     response   ...feel like harming myself
qu10        byte   %16.0g     response   ...feel like crying more than usual
qu11        byte   %16.0g     response   ...become annoyed or irritated easily
qu12        byte   %16.0g     response   ...have lost interest in other people
qu13        byte   %16.0g     qu13_t1    ...have trouble making decisions
qu14        byte   %16.0g     qu14_t1    ...feel unattractive
qu15        byte   %16.0g     qu15_t1    ...feel like not working
qu16        byte   %16.0g     qu16_t1    ...have trouble sleeping
qu17        byte   %16.0g     qu17_t1    ...feel tired or fatigued
qu18        byte   %16.0g     qu18_t1    ...makes my appetite lower than usual
qu19        byte   %16.0g     qu19_t1    ...concerned about my health
qu20        byte   %16.0g     qu20_t1    ...experience decreased libido

The responses consist of a 5-point Likert scale ranging from 1 (Strongly Disagree) to 5 (Strongly Agree). Questions 1-10 form the affective scale of the inventory and questions 11-20 form the physical scale. Data were simulated for 1000 imaginary people and included demographic variables such as age, sex and race. The responses can be summarized succinctly in a matrix of bar graphs:

Classical statistical analysis

The beginning of a classical statistical analysis of these data might consist of summing the responses for questions 1-10 and referring to them as the “Affective Depression Score” and summing questions 11-20 and referring to them as the “Physical Depression Score”.

egen Affective = rowtotal(qu1-qu10)
label var Affective "Affective Depression Score"
egen physical = rowtotal(qu11-qu20)
label var physical "Physical Depression Score"

We could be more sophisticated and use principal components to create the affective and physical depression score:

pca qu1-qu20, components(2)
predict Affective Physical
label var Affective "Affective Depression Score"
label var Physical "Physical Depression Score"

We could then ask questions such as “Are there differences in affective and physical depression scores by sex?” and test these hypotheses using multivariate statistics such as Hotelling’s T-squared statistic. The problem with this analysis strategy is that it treats the depression scores as though they were measured without error and can lead to inaccurate p-values for our test statistics.

Structural equation modeling

Structural equation modeling (SEM) is an ideal way to analyze data where the outcome of interest is a scale or scales derived from a set of measured variables. The affective and physical scores are treated as latent variables in the model resulting in accurate p-values and, best of all….these models are very easy to fit using Stata! We begin by selecting the SEM builder from the Statistics menu:

In the SEM builder, we can select the “Add Measurement Component” icon:

which will open the following dialog box:

In the box labeled “Latent Variable Name” we can type “Affective” (red arrow below) and we can select the variables qu1-qu10 in the “Measured variables” box (blue arrow below).

When we click “OK”, the affective measurement component appears in the builder:

We can repeat this process to create a measurement component for our physical depression scale (images not shown). We can also allow for covariance/correlation between our affective and physical depression scales using the “Add Covariance” icon on the toolbar (red arrow below).

I’ll omit the intermediate steps to build the full model shown below but it’s easy to use the “Add Observed Variable” and “Add Path” icons to create the full model:

Now we’re ready to estimate the parameters for our model. To do this, we click the “Estimate” icon on the toolbar (duh!):

And the flowing dialog box appears:

Let’s ignore the estimation options for now and use the default settings. Click “OK” and the parameter estimates will appear in the diagram:

Some of the parameter estimates are difficult to read in this form but it is easy to rearrange the placement and formatting of the estimates to make them easier to read.

If we look at Stata’s output window and scroll up, you’ll notice that the SEM Builder automatically generated the command for our model:

sem (Affective -> qu1) (Affective -> qu2) (Affective -> qu3)
(Affective -> qu4) (Affective -> qu5) (Affective -> qu6)
(Affective -> qu7) (Affective -> qu8) (Affective -> qu9)
(Affective -> qu10) (Physical -> qu11) (Physical -> qu12)
(Physical -> qu13) (Physical -> qu14) (Physical -> qu15)
(Physical -> qu16) (Physical -> qu17) (Physical -> qu18)
(Physical -> qu19) (Physical -> qu20) (sex -> Affective)
(sex -> Physical), latent(Affective Physical) cov(e.Physical*e.Affective)

We can gather terms and abbreviate some things to make the command much easier to read:

sem (Affective -> qu1-qu10) ///
(Physical -> qu11-qu20) ///
(sex -> Affective Physical) ///
, latent(Affective Physical ) ///
cov( e.Physical*e.Affective)

We could then calculate a Wald statistic to test the null hypothesis that there is no association between sex and our affective and physical depression scales.

test sex

( 1)  [Affective]sex = 0
( 2)  [Physical]sex = 0

chi2(  2) =    2.51
Prob > chi2 =    0.2854`

Final thoughts
This is an admittedly oversimplified example – we haven’t considered the fit of the model or considered any alternative models. We have only included one dichotomous independent variable. We might prefer to use a likelihood ratio test or a score test. Those are all very important issues and should not be ignored in a proper data analysis. But my goal was to demonstrate how easy it is to use Stata’s SEM builder to model data such as those arising from the Beck Depression Inventory. Incidentally, if these data were collected using a complex survey design, it would not be difficult to incorporate the sampling structure and sample weights into the analysis. Missing data can be handled easily as well using Full Information Maximum Likelihood (FIML) but those are topics for another day.

If you would like view the slides from my talk, download the data used in this example or view a video demonstration of Stata’s SEM builder using these data, please use the links below. For the dataset, you can also type use followed by the URL for the data to load it directly into Stata.

References

Beck AT, Ward CH, Mendelson M, Mock J, Erbaugh J (June 1961). An inventory for measuring depression. Arch. Gen. Psychiatry 4 (6): 561–71.

Beck AT, Ward C, Mendelson M (1961). Beck Depression Inventory (BDI). Arch Gen Psychiatry 4 (6): 561–571

Beck AT, Steer RA, Ball R, Ranieri W (December 1996). Comparison of Beck Depression Inventories -IA and -II in psychiatric outpatients. Journal of Personality Assessment 67 (3): 588–97
Bollen, KA. (1989). Structural Equations With Latent Variables. New York, NY: John Wiley and Sons

Kline, RB (2011). Principles and Practice of Structural Equation Modeling. New York, NY: Guilford Press

Raykov, T & Marcoulides, GA (2006). A First Course in Structural Equation Modeling. Mahwah, NJ: Lawrence Erlbaum

Schumacker, RE & Lomax, RG (2012) A Beginner’s Guide to Structural Equation Modeling, 3rd Ed. New York, NY: Routledge

Categories: Statistics Tags:

Stata YouTube channel announced!

StataCorp now provides free tutorial videos on StataCorp’s YouTube channel,

There are 24 videos providing 1 hour 51 minutes of instructional entertainment:

Stata Quick Tour (5:47)
Stata Quick Help (2:47)
Stata PDF Documentation (6:37)

Stata SEM Builder (8:09)
Stata One-way ANOVA (5:15)
Stata Two-way ANOVA (5:57)

Stata Box Plots (4:04)
Stata Basic Scatterplots (5:19)
Stata Bar Graphs (4:15)
Stata Histograms (4:50)
Stata Pie Charts (5:32)

And more are forthcoming.

The inside story

Alright, that’s the official announcement.

Last Friday, 21 September 2012, was an exciting day here at StataCorp. After a couple of years of “wouldn’t it be cool if”, and a couple of months of “we’re almost there”, Stata’s YouTube channel was finally ready for prime time.

Stata’s YouTube Channel was the brainchild of Karen Strope, StataCorp’s Director of Marketing, but I had something to do with it, too. Well, maybe more than something, but I’m a modest guy. Anyway, I thought it sounded like fun and recorded a few prototype videos. Annette Fett, StataCorp’s Graphic Designer, added the cool splash-screen and after a few experiments, we soon had 24 Blu-ray resolution videos. We’ve kicked off with videos covering topics such as a tour of Stata’s interface, how to create basic graphs, how to conduct many common statistical analyses, and more.

My personal favorite is the video entitled Combining Crosstabs and Descriptives because it’s relevant to nearly all Stata users and works well as a video demonstration.

Videos about Stata – isn’t that like dancing about architecture?

Stata has over 9,000 pages of documentation included in PDF format, a built-in Help system, and a collection of books on general and special topics published by Stata Press, and an extensive collection of dialog boxes that make even the most complex graphs and analyses easy to perform.

So aren’t the videos, ahh, unnecessary?

The problem is, it’s cumbersome to describe how to use all of Stata’s features, especially dialog boxes, in a manual, even when you have 9,000 pages, and 9,000 pages tries even the most dedicated user’s patience.

In a 3-7 minutes video, we can show you how to create complicated graphs or a sophisticated structural equation model.

We have three audiences in mind.

1. Videos for non-Stata users, whom we call future Stata users; videos intended to provide a loosely guided tour of Stata’s features.
2. Videos for new Stata users, such as the person who might simply want to know “How do I calculate a twoway ANOVA in Stata?” or “How do I create a Pie Chart?”. These videos will get them up and running quickly and painlessly.
3. Videos for experienced Stata users who want to learn new tips and tricks.

There’s actually a fourth group that’s of interest, too; experienced Stata users teaching statistics or data analysis classes, who don’t want to spend valuable class time showing their students how to use Stata. They can refer their students to the relevant videos as homework and thus free class time for the teaching of statistics.


text:
5 answers

# The Pade approximation of degree 5 is better than Taylor polynomial of degree 5.Select one:TrueFalse

###### Question:

The Pade approximation of degree 5 is better than Taylor polynomial of degree 5. Select one: True False

## Answers

#### Similar Solved Questions

1 answer
##### Case 6-30 Service Organization; Segment Reporting (L06-4) Music Teachers, Inc., is an educational association for music...
Case 6-30 Service Organization; Segment Reporting (L06-4) Music Teachers, Inc., is an educational association for music teachers that has 20,100 members. The association operates from a central headquarters but has local membership chapters throughout the United States. Monthly meetings are held by ...
1 answer
##### Design a 2-minute counter using JK Flip-Flops with every second equivalent to one clock cycle. Preferably...
Design a 2-minute counter using JK Flip-Flops with every second equivalent to one clock cycle. Preferably using Logisim....
5 answers
##### Wnal do Ihe entries of AQ represent? Select all that apply-DA The iotal cost of materials for economy models and deluxe models In Ihousands of dollars 08: The Iolal cost of labor and malerials (0t ( deluxo models in thousands of dollars Dc The Iolal cost ol Iahor lor economy modcls and deluxe modols thousands of dollars 0u The total cost ol Iabor and malorlal: Ior economy models In Ihousnd ol dollars
Wnal do Ihe entries of AQ represent? Select all that apply- DA The iotal cost of materials for economy models and deluxe models In Ihousands of dollars 08: The Iolal cost of labor and malerials (0t ( deluxo models in thousands of dollars Dc The Iolal cost ol Iahor lor economy modcls and deluxe modo...
5 answers
##### Find the inverse, f-1 of ff-1(x)= (Simplify your answer Use integers or fractions for any numbers in the The range of f(x) is (Type your answer in interval notation . Use integers or fractions for The domain of f-1 any (x) is (Type your answer in interval notation Use integers or fractions The range of f-1 for any (x) is (Type your answer in interval notation_ Use integers or fractions for any nu Enter your answer in the answer box and then click All parts Check Answer showing
Find the inverse, f-1 of f f-1(x)= (Simplify your answer Use integers or fractions for any numbers in the The range of f(x) is (Type your answer in interval notation . Use integers or fractions for The domain of f-1 any (x) is (Type your answer in interval notation Use integers or fractions The rang...
5 answers
##### Since X-X' and y-Y follows:@ @(x,y) = ar*yfur 0 < x <d,0 <y < dEnergy per unit length for Element 1:d2 + h2 Wen = zdh ( &1?+ a27 + Gagh+ aza3d+a3-We1 zdh 02 + 02 + 0 + 0 + 04) Wel 24? 0+4) 302
since X-X' and y-Y follows: @ @(x,y) = ar*y fur 0 < x <d,0 <y < d Energy per unit length for Element 1: d2 + h2 Wen = zdh ( &1?+ a27 + Gagh+ aza3d+a3- We1 zdh 02 + 02 + 0 + 0 + 04) Wel 24? 0+4) 302...
5 answers
##### Question 13; What would be the major product of the reaction below? (1R.2R)-1-ethoxy-_ methylcyclohexane EtOH 8. (1S 2R)-1-ethoxy-2-methylcyclohexane1-methylcyclohex- ene1-ethoxy-1-methylcyclohexaneQuestion 14: Is the reaction correct or not? Why or why not?~OkHzO , HzSOaReaction correct because the E2-like step is not reversibleB. Reaction is correct because of steric hindrarice by the cyclohexane. C. Carboxylic acid should've been fromed because the chromate ester is good leaving group_ D
Question 13; What would be the major product of the reaction below? (1R.2R)-1-ethoxy-_ methylcyclohexane EtOH 8. (1S 2R)-1-ethoxy-2-methylcyclohexane 1-methylcyclohex- ene 1-ethoxy-1-methylcyclohexane Question 14: Is the reaction correct or not? Why or why not? ~Ok HzO , HzSOa Reaction correct becau...
5 answers
##### The following scatterplot shows the relationship between anxiety test scores and science achievement test scores Choose the best interpretation of the relationship between anxiety level and science achievement:1 IAnxiety LevelThis graph shows strong negative linear relationship between anxiety and achievement in sdence_ This graph shows strong positive linear relationship between anxiety and achievement in sdence, This graph shows = moderate linear relationship between anxiety and achievement
The following scatterplot shows the relationship between anxiety test scores and science achievement test scores Choose the best interpretation of the relationship between anxiety level and science achievement: 1 I Anxiety Level This graph shows strong negative linear relationship between anxiety...
1 answer
##### (6 points) 10. Find the volume of the solid whose base is the region in the...
(6 points) 10. Find the volume of the solid whose base is the region in the xy-plane that is bounded by the parabola y = 4 – x2 and the line y=3x, while the top of the solid is bounded by the plane z = x + 6....
5 answers
##### MLtaFLa (LEZEne{Ne-tindDmn 0147KAernnnelHeLu4u4
MLta FLa (LEZEne {Ne-tind Dmn 0147 KAernnnelHe Lu4u4...
5 answers
##### If two resistors withresistances R1 and R2 areconnected in parallel, as in the figure below, then the totalresistance R, measured in ohms (Î©), is given by1R = 1R1 + 1R2. If R1 and R2 areincreasing at rates of 0.3 Î©/s and 0.2 Î©/s, respectively,how fast is R changingwhen R1 = 60 Î© and R2 = 90 Î©? (Roundyour answer to three decimal places.)
If two resistors with resistances R1 and R2 are connected in parallel, as in the figure below, then the total resistance R, measured in ohms (Î©), is given by 1 R = 1 R1 + 1 R2 . If R1 and R2 are increasing at rates of 0.3 Î©/s and 0.2 Î©/s, respectively, how fast is R changing when R1 = ...
5 answers
##### 3) Find the solution of the initial boundary value problemdtu = 02u, x â‚¬ (0,L), t > 0, u(z,0) = cos? (Tw/L); 81u(0,+) = &.U(L,t) = 0.Verify your answer. HINT Write cos?in terms of cos(z0) , cos(20)
3) Find the solution of the initial boundary value problem dtu = 02u, x â‚¬ (0,L), t > 0, u(z,0) = cos? (Tw/L); 81u(0,+) = &.U(L,t) = 0. Verify your answer. HINT Write cos? in terms of cos(z0) , cos(20)...
1 answer
1 answer
##### Please answer all parts and show work. 20) The equilibrium constant for e reaction Ni2+ (aa)...
please answer all parts and show work. 20) The equilibrium constant for e reaction Ni2+ (aa) H3 (ag) a Ni(NH3)62+ (ng) is th K-5.6x10s at 298 K. a) What is AG0 at this temperature?...
5 answers
##### Find Ihe Indicated score. Il graph right depicis IQ scores ults 100 and standard Ceviation ol 15 Ciick lo vler" pana eUlha_labla Clickto view pagu 2 oLuiha tablehose score:nontjIe disIrbulatTnriniStandard Normal Table (PogeTho indcPOSITIVE z Scores Stondjr Nortnal (DEubutlon' Keomt]E040Enlef you
Find Ihe Indicated score. Il graph right depicis IQ scores ults 100 and standard Ceviation ol 15 Ciick lo vler" pana eUlha_labla Clickto view pagu 2 oLuiha table hose score: nontjIe disIrbulat Tnrini Standard Normal Table (Poge Tho indc POSITIVE z Scores Stondjr Nortnal (DEubutlon' Keomt]...
5 answers
##### Wapplivd |e ubeti Me nebpcte bebrrmEatrardalkdteta binaltto undtns2n ( elilwgdAccording Ihe emr pirical tule . what percenlage 0l these (unes expecled 60 wribin +3 slandard deviabons the nean?Accordinqthe Chebyshev mule hat percentage of thesa funds are expactedt0 be Kthin 23 gtandard devialions In& Meant(Rounddecimal places needed )According " the Chebysher mule_aileast 93 7580 of these funds are expected - have cne-vear tatal relums betueen whattvo amounts?Between Irouna Io decimal
Wapplivd |e ubeti Me nebpcte bebrrm Eatrardalkdteta binaltto undtns2n ( elilwgd According Ihe emr pirical tule . what percenlage 0l these (unes expecled 60 wribin +3 slandard deviabons the nean? Accordinq the Chebyshev mule hat percentage of thesa funds are expactedt0 be Kthin 23 gtandard devialio...
5 answers
##### Find the sum of the series: k2 + k k=2a)2b)2d)3e)3
Find the sum of the series: k2 + k k=2 a) 2 b) 2 d) 3 e) 3...
1 answer
##### Your answer is partially correct. Try again. The following stockholders' equity accounts, arranged alphabetically, are in...
Your answer is partially correct. Try again. The following stockholders' equity accounts, arranged alphabetically, are in the ledger of Bridgeport Corporation at December 31, 2020. Common Stock (\$6 stated value) Paid-in Capital in Excess of Par-Preferred Stock Paid-in Capital in Excess of Stated...
1 answer
##### Please make sure that answer are available based on blue boxes below the problem How many...
Please make sure that answer are available based on blue boxes below the problem How many grams of NH Br must be dissolved in 1.00 L of water to produce a solution with pH = 5.16? The Kb of NH, is equal to 1.8 * 10- < PREV Based on your ICE table and Ka expression, determine the grams of NH,Br in...
1 answer
##### 4. Estimate the transition frequency for the poryphyrin molecule from m-11 to m 12, assuming that the pi electrons can be modeled as a particle in a ring of radius 440 picometers. (C 7. The...
4. Estimate the transition frequency for the poryphyrin molecule from m-11 to m 12, assuming that the pi electrons can be modeled as a particle in a ring of radius 440 picometers. (C 7. The most probable distance of the electron from the nucleus in a 1ls state hydrogen atom (with wavefunction V1) ca...
1 answer
##### (1Opoints) Find local extremne values &nd saclclle points of f(E,y) = +rv+1-34
(1Opoints) Find local extremne values &nd saclclle points of f(E,y) = +rv+1-34...
1 answer
##### 217 The space is C [0,21] with the inner product (f, g) = S fog(t)dt. Find...
217 The space is C [0,21] with the inner product (f, g) = S fog(t)dt. Find the third-order Fourier approximation to f(t) = 2 - 6. 0 The third-order Fourier approximation of fis-6 + 21 – 3 sint - 2 sin 2t - sin 3t| (Type an exact answer, using a as needed.)...
1 answer
##### Express as a sum or difference. $\sin (-4 x) \cos 8 x$
Express as a sum or difference. $\sin (-4 x) \cos 8 x$...
5 answers
##### C3 Thi s question Mhxoeult inletien7ice Light from [asct ic idcn (wo Vcnv nafTW' slls{EAAlu >JclJenenFuinc ( un the *crecil diteelly oposite the midpoint of the slits. (i) On 1*C; bclou- skctch thc variation with anglc â‚¬ of thc intensity ofthe light onthe scrccnJhicuaI(iip The scparation ofthc slits is 0.120 mm and thz wavelcugth ofthe lighc i~ 6.80 DI0' distance bctwccn Inc elits and thc screcn /> 4m Calculate the scparation 0f'thc bright fringss Uhe screenJ Sl A Cenedeik
C3 Thi s question Mhxoeult inletien7ice Light from [asct ic idcn (wo Vcnv nafTW' slls {EAAlu >Jcl Jenen Fuinc ( un the *crecil diteelly oposite the midpoint of the slits. (i) On 1*C; bclou- skctch thc variation with anglc â‚¬ of thc intensity ofthe light onthe scrccn JhicuaI (iip The scp...
5 answers
##### I1 L 1 1 1 8 1 0 1 ! 7 8 LL 1 d 1 1 193 1 6 21 K 1 2 1 2 1 4 1 ! 1 0 4 V L L 1 Ww li 7 4 [ ] 4 1 Ub 9 4 1 1 4 1 2 1 1 ] 2 1 4 1 1 1 1 1 W 4 3 4 9 I 1 1 1 I 1! 3 0] 1 1 9{ 1 J1 3
I1 L 1 1 1 8 1 0 1 ! 7 8 LL 1 d 1 1 193 1 6 21 K 1 2 1 2 1 4 1 ! 1 0 4 V L L 1 Ww li 7 4 [ ] 4 1 Ub 9 4 1 1 4 1 2 1 1 ] 2 1 4 1 1 1 1 1 W 4 3 4 9 I 1 1 1 I 1! 3 0] 1 1 9 { 1 J1 3...
5 answers
##### Complete the table by writing the alkene and reaction conditions required: ALKENE REACTANTS/CONDITIONS PRODUCT0-H
Complete the table by writing the alkene and reaction conditions required: ALKENE REACTANTS/CONDITIONS PRODUCT 0-H...
5 answers
##### ContenSelect the During Describe [ Cnntennattibution any liquid Bas solid Question of the correct homogeneous homogeneous above anstei 1 below: catalysis; catal 4 present in the liquid phase; FEEDBACK which phase MORE must the 1 reactants1 37
conten Select the During Describe [ Cnntennattibution any liquid Bas solid Question of the correct homogeneous homogeneous above anstei 1 below: catalysis; catal 4 present in the liquid phase; FEEDBACK which phase MORE must the 1 reactants 1 3 7...
1 answer
##### Discuss in detail three factors which motivated European contact, conquest, and colonization of non-western societies in...
Discuss in detail three factors which motivated European contact, conquest, and colonization of non-western societies in Africa, Asia, and Latin America....
1 answer
##### 5. Inflation and CPI Suppose in Economy A, the final goods produced in year 2015, 2016...
5. Inflation and CPI Suppose in Economy A, the final goods produced in year 2015, 2016 and 2017 are milk, beef, gasoline, cloth. The quantities (Q) produced and prices (P) are given by the following: Year 2015 Qmilk Pmilk 200 2.35 220 2.49 233 2.79 (beer Pbeer lgasolin 120 3.30 400 123 3.35 430 145 ...
1 answer
##### STRATEGIC MANAGENMENT
Question 1  (a) Assume that the spot exchange rate of the South African Rand is 1.73.   How will this spot rate adjust according to Purchasing Power Parity (PPP) if South Africa experiences an inflation rate of 7 percent while Zambia experiences an inflation rate of 2 percent? &n...
5 answers
##### Y3 4 22y5 _ eUse implicit differentiation t0 find formula forsin (&)2e2 x "2ws 3 y2 xSy5Your answer can be function of both and y:HintPenaltyNotation0.0view AgainWarning 0.0viaw Again
Y3 4 22y5 _ e Use implicit differentiation t0 find formula for sin (&) 2e2 x "2ws 3 y2 #xSy5 Your answer can be function of both and y: Hint Penalty Notation 0.0 view Again Warning 0.0 viaw Again...

-- 0.218094--


text:
# Logistic Regression Model in R

Logistic regression aims to model the probability of an event occurring depending on the values of independent variables.

The logistic regression model seeks to estimate that an event (default) will occur for a randomly selected observation versus the probability that the event does not occur. Suppose we have data for 1000 loans along with all the predictor variables and also whether the borrower defaulted on it or not. Here the probability of default is referred to as the response variable or the dependent variable. The default itself is a binary variable, that is, its value will be either 0 or 1 (0 is no default, and 1 is default).

In logistic regression, the dependent variable is binary, i.e. it only contains data marked as 1 (Default) or 0 (No default).

We can say that logistic regression is a classification algorithm used to predict a binary outcome (1 / 0, Default / No Default) given a set of independent variables. It is a special case of linear regression when the outcome variable is categorical. It predicts the probability of occurrence of a default by fitting data to a logit function.

R is a versatile package and there are many packages that we can use to perform logistic regression. In this case study we will use the glm() function in R. R also has a very useful package called caret (short for classification and regression training) which streamlines the process of training models using different algorithms.) We will use the caret package in the second case study.

We will train the model using the training data set, credit_train. Just to keep the whole process simple, we will use only a handful of selected variables in the model. To start with let's just use five variables to determine the value of Creditability.

> set.seed(1)
> LogisticModel <- glm(Creditability ~ Account.Balance + Payment.Status.of.Previous.Credit + Purpose + Length.of.current.employment + Sex...Marital.Status, family = binomial, data = credit_train)

Let’s look into what this model contains:

> LogisticModel
Call:  glm(formula = Creditability ~ Account.Balance + Payment.Status.of.Previous.Credit +
Purpose + Length.of.current.employment + Sex...Marital.Status,
family = binomial, data = credit_train)
Coefficients:
(Intercept)                    Account.Balance2
-2.1396                              0.5810
Account.Balance3                    Account.Balance4
1.0231                              1.8129
Payment.Status.of.Previous.Credit1  Payment.Status.of.Previous.Credit2
0.2335                              1.1462
Payment.Status.of.Previous.Credit3  Payment.Status.of.Previous.Credit4
0.9977                              1.8753
Purpose1                            Purpose2
0.8587                              0.3922
Purpose3                            Purpose4
0.6738                              0.2800
Purpose5                            Purpose6
0.4248                             -0.5447
Purpose8                            Purpose9
1.9817                              0.1309
Purpose10       Length.of.current.employment2
0.3163                             -0.1731
Length.of.current.employment3       Length.of.current.employment4
0.1855                              0.5960
Length.of.current.employment5               Sex...Marital.Status2
0.2057                              0.2572
Sex...Marital.Status3               Sex...Marital.Status4
0.4461                              0.5433
Degrees of Freedom: 699 Total (i.e. Null);  676 Residual
Null Deviance:        871.5
Residual Deviance: 723.3     AIC: 771.3
>

Mathematically this can be represented as follows:

In our equation, B0 = -2.1396, B1 = 0.5810, and so on.

We can now use this model to predict the value of Creditability for new data. Let’s say we have a new applicant. We can predict the Creditability of this applicant. Let’s create some random features data for this new applicant:

P(Creditability = 1 | Account.Balance = 4, Payment.Status.of.Previous.Credit = 2, Purpose = 1, Length.of.current.employment = 4,Sex…Marital.Status = 2)

\= 1 / ( 1 + e^-( -2.1396 + 1.8129*1 + 1.1462*1+ 0.8587*1+0.5960*1+0.2572*1))

\= 0.9263139694

We can also use the predict() function to do the same job. Formula: predict(model, new.df)

> newapplicant <- data.frame(Account.Balance=as.factor(4), Payment.Status.of.Previous.Credit=as.factor(2), Purpose=as.factor(1), Length.of.current.employment=as.factor(4), Sex...Marital.Status=as.factor(2))
> result <- predict(LogisticModel, type = 'response', newdata = newapplicant)
> result
1
0.9263105
>

Response gives you the numerical result while class gives you the label assigned to that value. Response lets you to determine your threshold. We can determine a threshold value say 0.6, and say if result > 0.6, credibility = 1, else credibility = 0.

> if(result>0.6) {credibility = 1} else {credibility = 0}
> credibility
[1] 1

That done, we can move on to fitting the model we just created to the test set, credit_test, and prepare to make our first prediction.

predicted_values <- predict(LogisticModel, type = 'response', newdata = credit_test)

We have successfully created and fitted our first model. The predicted_values is a numeric list of the predicted value for each borrower in the credit_test dataset.

> plot(predicted_values)

### Class Labels

The predicted values are in the range of 0 and 1. We can apply the cutoff (-.6) and apply labels (0 and 1) to all the predicted values.

> pred_value_labels = rep(as.factor("0"),length(credit_test))

> pred_value_labels = rep("0",length(credit_test[,1]))
> pred_value_labels[predicted_values>.6] = "1"
> pred_value_labels <- as.factor(pred_value_labels)
> pred_value_labels
[1] 1 0 1 1 1 1 1 0 1 1 1 0 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 0 0 0 0 1 0 1 1 1 0 0 1 1 0 0 0 1
[50] 0 0 1 0 1 1 1 0 1 1 0 0 0 0 1 0 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1
[99] 1 0 1 1 0 1 1 1 1 0 1 0 1 1 1 0 0 0 0 1 1 1 1 0 0 0 1 1 0 0 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 0 1 1 0
[148] 0 1 1 1 1 0 0 1 0 1 1 1 0 1 0 1 1 1 0 1 0 0 0 0 0 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 0 1 1 1 1 1 1 1 0
[197] 1 1 0 1 0 1 1 1 0 1 0 1 0 0 0 1 1 0 1 0 0 0 1 1 1 1 0 1 1 0 0 1 1 1 0 0 1 1 1 1 1 1 1 1 0 0 1 1 0
[246] 1 1 0 1 1 0 0 1 1 1 1 1 0 0 1 1 1 0 0 0 1 1 1 0 0 1 0 1 0 1 1 1 0 1 1 0 1 0 1 1 1 1 0 1 0 0 1 0 0
[295] 1 1 1 0 0 0
Levels: 0 1
>

In the next lesson we will learn about how we can statistically measure the performance of the fitted model.


text:
View Single Post
 Emeritus Sci Advisor Thanks PF Gold P: 15,868 And for discontinuous functions, the result might not even be true!! It needs to be: If $\int_A f(x)dx = \int_A g(x)dx$ then $f=g$ almost everywhere.