title: Gradient Based Learning Algorithms vs Global Optimization Learning Algorithms for Neural Networks

text: Neural Networks are usually trained using a gradient based learning algorithm, such as the back propagation algorithm or some variant of it, but can you use global optimization algorithms, such as the Genetic Algorithm, Nelder-Mead Polytope Algorithm, and Particle Swarm Optimisation to train the network? Since training a neural network all boils down to minimizing a multi-variable cost function, I would assume that it should be easy to do this using a global optimization method, but I have tried to do it myself and I'm getting very bad results. The GA algorithm reduces my cost function from around 250 (when it is input with random synaptic weights), to only around 170. The Nelder-Mead Algorithm would apparently take years, and I have not yet tried PSO as there is no inbuilt MATLAB function for it. So is it just accepted that gradient based algorithms are the most suitable for training a Neural Network? If so can someone please point me towards a source of this information? This will be very helpful as I could reference the source in my project to explain why I gave up trying to use global optimization methods to train the network. Thanks!


title: Newton's algorithm for a Minimax Procedure

text: In order to find the minimax solution of a certain problem I have to determine the constant $c$ that provides the root for the following equation: $$l(c)=3\times \left [ 1-\Phi \left( c-75 \right) \right]-\Phi \left(c-78 \right)=0$$ where $\Phi (.) $ is as usual the CDF of a standard normal random variable. Using Newton's algorithm then we have to iterate the following: $$c_1=c_0-\frac{l(c_0)}{l \prime (c_0)}$$ where $l \prime (c)=-3 \times \phi \left(c-75 \right)-\phi \left(c-78 \right)$ and the corresponding R-code for say, 100 iterations: > for(i in 1:100){ > > l=3* (1-pnorm(c-75))-pnorm(c-78) > > lprime=-3* dnorm(c-75)-dnorm(c-78) > > c=c-(l/lprime) > > } > > c I can see from the book that the solution should be $76.8$. The problem is that the algorithm reaches that value only for a small neighbourhood around $c=75$ and does not converge otherwise (try it for $c=60$ or $c=80$). The algorithm is then essentially useless as in most cases I would not have knowledge of the solution and therefore I would have to pick a very rough initial estimate. Is there anything I can do to improve the algorithm? If not, is there perhaps another algorithm that converges for a larger set of initial values? I know this is more of a mathematical question but I would be very interested in hearing the opinions of people with experience in statistics and numerical optimization so please do not move it, at least not right away. Thank you.


title: Percentile Loss Functions

text: The solution to the problem: $$ \min_{m} \; E[|m-X|] $$ is well known to be the median of $X$, but what does the loss function look like for other percentiles? Ex: the 25th percentile of X is the solution to: $$ \min_{m} \; E[ L(m,X) ] $$ What is $L$ in this case?


title: Validity of normality assumption in the case of multiple independent data sets with small sample size

text: Due to limitations in experimental setup, I only have small data sets with n=3. Despite the low df the difference between treated and control is large enough to generate a significant p-value. The problem is that with small sample sizes doing a t-test becomes more sensitive to the assumption that the data are drawn from a population of a normal distribution. In my case however, multiple independent experiments consistently yield a similar result. I cannot group the data of the experiments because of small differences in the data. For example the increase between treated and control in one experiment is slightly bigger than in another experiment, which is likely caused by small variances in experimental conditions (it would get technical to explain this further). Despite this the same increase is consistently observed and each time the 3 data points have a small standard deviation for both groups. So my question is whether it is defensible to make the normality assumption based on the data of multiple independent experiments with a small sample size? If not am I right that it would not be appropriate to use any statistics in this case?


title: Bias of maximum likelihood estimators for logistic regression

text: I would like to understand a couple of fact on maximum likelihood estimators (MLEs) for logistic regressions.   1. Is it true that, in general, the MLE for logistic regression is biased? I would say "yes". I know, for example, that sample dimension is related to the asymptotic bias of MLEs.  Do you know any elementary examples of this phenomenon?   2. If the MLE is biased, is it true that the covariance matrix of the MLEs is the inverse of the Hessian of the maximum likelihood function?  _edit_ : I have met this formula quite often and without any proof; it seems a quite arbitrary choice to me.


title: IQR for Normal Probability Distribution

text: I was attempting to find the IQR for a normal distribution with mean being 4.5, while the standard deviation. My attempt yielded the result **2.144** , while the actual answer provided was **2.1584**. Following is my workings, appreicate if anyone will be able to guide me where I am going wrong please: > X ~ N(4.5, 1.6^2) > > P(A < X < B) = 0.5 > > thus using table 4 allows us to deduce that Za = -0.67, while Zb = 0.67. > > using the formulae x = z . Population sd + population mean, > > A = (-0.67)(1.6) + 4.5 = 3.428 > > B = (0.67)(1.6) + 4.5 = 5.572 Thus, for IQR, I took 5.572 - 3.428 = 2.144.


title: Can (should?) regularization techniques be used in a random effects model?

text: By regularization techniques I'm referring to lasso, ridge regression, elastic net and the like. Consider a predictive model on health care data containing demographic and diagnosis data where length of stay for inpatient stays is being predicted. For some individuals there are multiple LOS observations (i.e., more than one IP episode) during the baseline time period which are correlated. Does it make sense to build, for example, an elastic net predictive model which contains a random effect intercept term for each individual?


title: Interpreting Labour Economics regression model

text: Instrumental variables (2SLS) regression   Number of obs = 603 F( 5, 597) = 41.96 Prob > F = 0.0000 R-squared = 0.1386 Root MSE = .26523                        |                     logwage |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]            higher |   .2566441   .1030163     2.49   0.013     .0543257    .4589625          age |  -.0054228   .0055142    -0.98   0.326    -.0162524    .0054068         age2 |   .0000223   .0000675     0.33   0.741    -.0001103    .0001549     urban    |   .2587201   .0304102     8.51   0.000     .1989962    .3184441        rural |  -.0207374   .0296745    -0.70   0.485    -.0790164    .0375417        _cons |   2.370127   .1065349    22.25   0.000     2.160899    2.579356      Instrumented: higher ## Instruments: age age2 metropolitan rural father_educated I am doing an assignment titled does expansion of higher education improve the earnings? The case of Russia Regression has been run as shown above, the regression model is based on the widely used Mincer Model with minor modification and proxy to it. Can anyone of you who are kind enough to assist in term of what do you think of the model and its output? And how would you interpret the output?


title: Minimum Distance Estimator in r

text: I am running a regression to estimate the relationship between a dependent variable and some regressors. Let's say I have two alternative models based on two different assumptions about the exogeneity of the regressors involved. Therefore I get two set of parameters, say `ß` and `b`, both of them including `n` parameters for the `n` regressors. I have some theoretical justification to say that the former set is more likely to be obtained from the correct specification of the model. However, I would like to provide a third set of parameter, say `a`, obtained as the Minimum Distance Estimator, according to the reasoning proposed for instance in Green MDE §13.3 . Is there any R package allowing me to obtain my third vector of parameters ? I found the actuar package that should be doing the job I need. I understand all the values to be specified in the function `mde`, still I am a bit confused about the value `fun`. What am I exactly doing by specifying the `fun` argument ? Have I to repeat the procedure for each parameter or is there a way to do that for all the `n` parameters I estimated ?


title: Finding structural breaks in heteroskedastic time series

text: I'm trying to identify structural breaks in the movement of reserve currencies. I'm not yet all that versed in the finer details of time series, but I've been reading up on ARCH and GARCH estimators. However, in my model I don't really want to look at any external factors, like government interventions, gdp, macroeconomic drivers, etc. Instead I want to establish a certain confidence interval (on returns or variance?) and identify abnormal movement outside that bandwidth. The classical structural break identifiers seem to me to be very "exclusive", meaning they are good at identifying massive underlying changes; however for me it wouldn't be a problem to have some false positives, because I want to test how breaks over a certain magnitude occur. I was thinking of employing some sort of rolling vintage, or segmenting the data at different rates so that changes can be identified, but my mind is all over the place in this regard. What I'm asking is how you would go about detecting abnormal variance or behavior in a univariate timeseries with tons of datapoints. Do I have to go and set up an ARCH/GARCH model or is there a simpler, cruder way for an algorithm that I can tune to give me points in which behavior of the time series changes?


title: How to extract ngrams from ambigous text after lemmatization?

text: After lemmatization of text I have a sequence of sets of lemmas, because every word can correspond to more than one lemma. How should I extract ngram statistics based on that? The only thing that comes to mind is counting every possible combination, but this seems a bit weird. Is there any other option? More implementation-related question: is there a way to do this by NLTK?


title: Propensity score in time-updated (counting process) Cox regression

text: Suppose I aim to examine if treatment A is better than treatment B. My dataset consists of 10,000 individuals with a total of 100,000 observations. Thus each individual is observed 10 times and his covariates (blood pressure, blood lipids etc) are updated each visit. To provide more robust estimates, I use all observations in the Cox regression, by counting process syntax. I would like to adjust for propensity score as well but the problem is: (1) Should the propensity score be calculated for each observation (i.e 10 times for each individual)? (2) Should the propensity score be calculated once (per individual) and then held fixed? **My second question** I have read R D'Agostinos review (http://www.stat.ubc.ca/~john/papers/DAgostinoSIM1998.pdf) on use of propensity score but after reading a paper (which investigated differences in treatment A and treatment B on risk of death), where authors stated that their covariance adjustment was done only by stratification with a propensity score; is it not "better" to include the covariates in the Cox regression as well, despite being included in the prop score? I did not attach a data set because this is a methodological discussion. Greetings


title: Backward selection in a Cox regression model

text: My goal is to fit a cox regression model in SAS, for which I use the `PROC PHREG` statement. As I am still new to regression methods, I would appreciate a little of your help. My procedure is as follows: First, I check all variables in a univariate regression and select those, which have a p-value less than 0.25. The next step is that these significant variables (p-value < 0.25) are considered in a model with backward selection. Variables with a p-value less than 0.167 (according to the AIC criterion) are removed from the model. The reduced model should then include the best explanatory variables. My question now is, whether I can find out, why one variable is removed in the backward selection or which variable rather explains the influence? I am surprised, that this variable drops out, although it is very significant in the univariate regression. Furthermore, I have deleted a variable from the backward selection, as this one is removed later anyway, and see, my "problem variable" remains in the model. How does this work? Thank you for your thoughts!


title: Spotting trends in time based data

text: I have a dataset which I need to spot trends in. The actual data refers to operation calls which take a certain amount of time to complete. My client wants to know which operation calls are improving and which are deteriorating over a time period. Each operation call has an entry for each day which shows the average time that call took to return.               Date        Call    Avg duration     -----------------------------------     18/03/2013  Call 1  123ms     19/03/2013  Call 1  156ms     20/03/2013  Call 1  198ms     21/03/2013  Call 1  99ms      My current implementation uses               (CurrentDayAvgDuration - PrevDayAvgDuration) / √(CurrentDayAvgDuration + PrevDayAvgDuration)      for each day and then takes the average of that value over all days but I'm not too sure if the value is relevant. That formula is one I found and plugged in just to try. My question is, how can I calculate a single value per operation call which can be used in comparison with all of the other operation calls to determine which are trending upwards (deteriorating) and which are trending downwards (improving)?


title: Frailty model using bootstrap sampling

text: Let say I have a 'kidney catheter' data set which is built in `survival` package. Data are about the recurrence times to infection, at the point of insertion of the catheter, for kidney patients using portable dialysis equipment. Catheters may be removed for reasons other than infection, in which case the observation is censored. Each patient has exactly 2 observations. **We can use shared frailty model to model this correlated survival times. The unobserved risk factor that is common for all the observations in the same individual or cluster is called the frailty. By performing the bootstrap, we obtain for each simulated data set a bootstrap estimate of the heterogeneity or frailty parameter. How can we estimate the frailty parameter with model based bootstrap approach using`R code`?**


title: Confidence intervals for extreme value distributions

text: I have wind data that i'm using to perform extreme value analysis (calculate return levels). I'm using _R_ with packages 'evd', 'extRemes' and 'ismev'. I'm fitting GEV, Gumbel and Weibull distributions, in order to estime the return levels (RL) for some period T. For the GEV and Gumbel cases, i can get RL's and Confidence Intervals using the _extRemes::return.level()_ function. But for the Weibull case, i need some help to generate the CI's. Any ideas? Thanks!


title: ANOVA in R: How MSresiduals are estimated when using Error(interaction)?

text: I have checked the existing answers, but I found no answer to my question. It is about choosing the right ANOVA model in R and how the MSresiduals are estimated. I have an experiment with two factors: `MANAGEMENT` and `REGION`. I have three `MANAGEMENT` types and three `REGIONS`. That means nine sites under study but no repetition, like this: `REGION`: `A`, `B` and `C` (those regions are not clonal copies of each other)   `MANAGEMENT`: `N`, `P`, `ST`   the dependent variable is called `PHt` I have set the `aov` in R in two possible ways. Model1: `aov(PHt~REGION+MANAGEMENT+REGION*MANAGEMENT, data=data)` Model2: `aov(PHt~REGION*MANAGEMENT+Error(subject/(REGION*MANAGEMENT)),data=data)` In respect to these two models I have three questions:   1. Is it correct that the second model should be used if we aim to include the interaction term into the `Error` term?    2. How does the model 2 estimate the MSresiduals? It seems to use different MSresiduals for each variable. Why? (See output for both models below.)   3. Since I have unbalanced data, is R using an approximation of a Wald test? Here are results of the two models.               Model 1     Model <- aov(PHt~REGION+MANAGEMENT+REGION*MANAGEMENT, data=data)     summary(model)                         Df      Sum Sq        Mean Sq    F value       Pr(>F)         REGION               2      0.05          0.0229       0.907        0.4039         MANAGEMENT           2      0.18          0.0884       3.494        0.0306 *       REGION:MANAGEMENT    4      3.34          0.8351      33.021        <2e-16 ***     Residuals           1744    44.10         0.0253                        Model 2     aov(PHt~REGION*MANAGEMENT+Error(subject/(REGION*MANAGEMENT)),data=data)          SUMMARY:       Df Sum Sq Mean Sq F value Pr(>F)     Residuals 195  4.714 0.02494                         Error: subject_f:REGION                Df Sum Sq Mean Sq F value Pr(>F)     REGION      2  0.066 0.03281   1.406  0.246     Residuals 388  8.822 0.02334                         Error: subject_f:MANAGEMENT                 Df Sum Sq Mean Sq F value Pr(>F)       MANAGEMENT   2  0.205 0.10260   4.413 0.0127 *     Residuals  388  8.788 0.02325                      Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Error: subject_f:REGION:MANAGEMENT                        Df Sum Sq Mean Sq F value Pr(>F)         REGION:MANAGEMENT   4  3.451  0.8627   31.99 <2e-16 ***     Residuals         776 20.389  0.0270                        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


title: example where comparison of three mean anova and t-test have different results

text: suppose to have sample from 3 groups A,B,C. The hypothesis `H0: the mean of the 3 groups is the same` can be tested using 3 independent t test.               test1: mean(A)=mean(B) level 0.05     test2: mean(B)=mean(C) level 0.05     test3: mean(A)=mean(C) level 0.05      It's known that we should prefer an ANOVA test because with the previous method there is an increasing risk of type 1 error. I would like to have an example with simulated data where the first procedure lead us to an error while the ANOVA return the good result. The best would be an R code to simulate the experiment. Thanks!


title: parametric bootstrap for low sample sizes

text: I believe that this question is sufficiently different from previous related ones to warrant a new post. (I apologize if it has been answered already) I need to decide between various resampling methods to "best" (highest power and correct type-I error to reject H0 for the right reasons) evaluate auto correlation in binary sequences of small to moderate lengths (20-50) with relatively low incidence rate. The time series are too short to justify any normal approximations, so to e.g. to decide whether or not the observed auto correlation is "real", I would like to generate my empirical Null distribution under the assumption of no auto correlation. With that goal in mind I can either   1. shuffle the sequence repeatedly (fixed number of 0s and 1s), or   2. estimate the probability $\hat{p}$ of an event and repeatedly draw from the respective Bernoulli process. The latter is -I think- known as the parametric bootstrap. I struggle with the correct choice for the following reasons: (i) The estimate for the true Bernoulli probability will be rather poor/noisy for low smple sizes. That bias will be part of the generated empirical NULL. (ii) The permutation procedure 1 will generate test statistics with considerably less variation than procedure 2. In fact, the distribution will tend to generate few discrete levels. What line of reasoning will defend either choice ? Here is an example:               x=c(0,1,0,0,0,1,1,0,0,1,0,0,0,0,1,0)     #observed acf:     aObs = acf(x,lag.max=1,plot=F)$acf[2]         a1=a2=rep(NA,100)         for (i in 1:100){           a1[i]=acf(sample(x),lag.max=1,plot=F)$acf[2]           a2[i]=acf(rbinom(length(x),1,mean(x)),lag.max=1,plot=F)$acf[2]     }     hist(a1);abline(v=aObs,lty=2,col=2)     hist(a2);abline(v=aObs,lty=2,col=2)


title: Kullback-Liebler Divergence function in R package "FNN" returns NAN for certain data? Why?

text: I am using the "FNN" package to calculate the Kullback-Liebler Divergence between two numeric vectors:               require(FNN)     load("~/xy.RData")     KL.divergence(x,y)      (Link to the data file is given at the bottom) The result is:               NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN      The two numeric vectors are clean (no missing data). The histogram of x is: ![enter image description here](http://i.stack.imgur.com/muIgq.png) And the histogram of y is: ![enter image description here](http://i.stack.imgur.com/IEcug.png) x contains 1280 entries, and y contains 7824 entries. For comparison:               KL.divergence(rnorm(1280),rnorm(7824))      works fine and yields real valued numbers. Both x & y contain a lot of zeros, 1s, and -1s. So I tried removing those:               KL.divergence(x[ (x!=1) & (x!=0) & (x!=-1) ],y[ (y!=1) & (y!=0) & (y!=-1) ])      But I still get all NaNs for results. So maybe something specific to the data set is causing this? Or maybe there is a bug with KL.divergence? Any insight would be appreciated! EDIT: If any one is interested, here is the data file: http://www.wikiupload.com/BZGLYFEGZRYQJXN (hit the green button that says "Download This File")


title: Use of the Jeffreys prior in multidimensional models

text: Suppose a model, $$x_{i} \sim N(\theta_{i}, \phi), \text{ for } i=1,\ldots,n$$ Furthermore, suppose the variance parameter, $\phi$, is some known constant. The multidimensional Jeffreys prior is defined as, $$\pi(\vec{\theta}) \propto \sqrt{\text{det}I(\vec{\theta})}$$ where $\text{det}I(\vec{\theta)}$ is the determinant of the multidimensional Fisher information. Can a multidimensional Jeffreys prior be used to formulate a prior distribution for the posterior distribution, $p(\vec{\theta}|x_{1},\ldots,x_{n})$?


title: Sample changing distribution

text: I have the following process: I have N buckets in front of me, $M$ of which are filled with water (the other ones are empty).   1. I pick one of them (consider uniform distribution) and empty it (so now I have $M-1$ of them filled with water).    2. Afterwards $P$ more get filled with water (picked by some unknown rule) (so I have $P+M-1 < N$ filled buckets).    3. go to 1 How could I model it to be as fair as possible with the picking of the bucket at each step? (Let's say we give some preference to buckets who had been filled earlier than others have.)


title: Three-prisoner problem and Bayes rule

text: Here is the wiki of Three-prisoner problem, in which only one prisoner is pardoned, and the Bayes solution is given in the wiki. My problem is pretty much the same, except that only one prisoner is executed. Here is the description: > A, B, C are three prisoners, two of which will be pardoned, A asks the > warden to tell him which ONE will be pardoned other than himself. A thinks > that before asking, he has a chance of $\frac{2}{3}$ to be pardoned, but if > he asks, the chance goes down to $\frac{1}{2}$. A's calculation is wrong, > explain why. I wanted to do it the Bayes way, let $A,B,C$ be the event that A,B,C will be pardoned repectively, and let $b$ be the event that the warden mentions B will be pardoned. The prior here is $P(A) = \frac{2}{3}$, and after observing $b$, the posterior of $A$ is $$P(A|b) = \frac{P(b|A)P(A)}{P(b)}$$ . I think $P(b|A)P(A) = \frac{1}{2}\frac{2}{3}$, since if $A$ happens, the warden should not mention $A$, he should mention the other one who will be pardoned, which has probability $\frac{1}{2}$ for B. And $$P(b) = P(b|A)P(A) + P(b|B)P(B) + P(b|C)P(C)$$ , right? How to calculate this?


title: Comparing two points corresponding to two different normal distributions

text: I have two multi variate normal distributions N1 and N2. Say two points p1 is from N1 and p2 is from N2. I want to get some statistical features from these two points. How can I do it? I need a richer representation (not one number). Say my points are in a d dimensional space, in that case I need a d dimensional vector which represents the statistical difference between these two points.


title: Why does SPSS throw out variables in OLS (and how can I prevent this from happening)?

text: My dataset consists of 4,000 companies from 24 countries. Each company has certain features, the IV's. Consider the variable list below: green is the DV, red the IV's and blue are dummies (one for each country a company might be from) for which I want to control the OLS output. There are 24 countries entered. ![enter image description here](http://i.stack.imgur.com/KLz3z.png) However, once I let SPSS perform an OLS on these variables, it throws out most of the IV's. ![enter image description here](http://i.stack.imgur.com/qlpWP.png) The resulting model is almost useless in explaining any of the variance in the DV. My questions are:   1. Why is SPSS throwing out these variables?   2. How can I prevent this from happening? ![enter image description here](http://i.stack.imgur.com/iwd0K.png) `[edit]` Added regression model in response to @Penguin_Night's and @NickCox's answers: ![enter image description here](http://i.stack.imgur.com/UxFK3.png)


title: Which is a good tool to compute parameters for a beta-distribution?

text: I have a β-distribution(which looks like a normal distribution), for which I want to calculate the values of α and β using maximum likelihood estimation method. I am new to statistics, so would like to know if there is any tool available which can be used to compute these parameters. PS: Similar question has been answered before, but the solutions were focused towards using R statistical package. I am looking for an out of the box solution, or a solution I can code in a standard scripting language like python.


title: p-value of t-test versus F-test(joint hypothesis)

text: SAS code :               data pr;         input y x1 x2 @@;         cards;         5 2 2 7 3 3 3 2 0         5 2 4 4 3 3 7 2 4         ;       run;       proc reg data = pr;         model y = x1 x2; run; quit;        When I run the regression procedure, I got the result               [F-test : p-value for model = 0.3757]       [t-test : p-value for x1 = 0.9238]       [t-test : p-value for x2 = 0.2045]        I'm wondering why the p-value of F-test is bigger than that of t-test for x2? F test has the joint hypothesis H0 : beta1 = beta2 = 0   so, it has the more complicated hypothesis than t test (H0 : beta2 = 0) My logic is this:   when hypothesis has the more restrictions, then it is far difficult to realize, so p-value gets lower. Where am i misunderstood?


title: Econometrics : Multiple regression Fisher and Student statistics

text: I am trying to estimate a production function called Cobb-Douglas. For the period 1958 to 1972 and for the agricultural sector in Taiwan, we observed:               An: Year of observation     Y: Real production in millions New Dollars Taiwan (NDT)     L: Day of Labour, in millions     K: Real Capital in millions de NDT      We consider the following model M1: `lm(formula = LY ~ LL + LK)` for = 1958, ..., 1972 where **the variables have been log-transformed.**   1. I try to test H0: $\beta_{1}+\beta_{2}=1$ versus H1: $\beta_{1}+\beta_{2}\neq 1$ with a Fisher statistic (bilateral test)   2. I try to test the same hypothesis with a Student statistic this time.   3. I try to test (unilateral test) H0: $\beta_{1}+\beta_{2}\leq 1$ versus H1: $\beta_{1}+\beta_{2}>1$ with a Student statistic


title: Unique time variable panel regression fixed effect

text: Currently I am working on my master thesis which is about risk-adjusted returns (Sharpe ratio) of Asian REITs. I just transformed all the data in variables which are ready to use in Stata. My data consist of year (2002 till 2012), Key (cusip identifier), Country, and all sorts of firm specific characteristics. For declaring a panel data in Stata, we have to declare two variables - Panel and Time. But Time variable must be unique, that is, no duplicate time in each panel. But when I declare "Country" as panel variable and "Year" as time variable, there is a problem. Problem is the duplicate year in each panel. For example, Japan has several entries for the year 2002. But for panel data there should not duplicate time values. I think I need to define a unique time variable or drop time variable from the analysis. Using only panel variable I can run panel data regression, but time series analysis such as unit root can't be run. So, question is how I declare a time variable?


title: Linear regression on several trials of experiment

text: Sorry about asking such a basic question. Assuming I have data like this `x Trial1 Trial2 Trial3 1 1.0 2.0 3.0 2 1.1 2.1 3.1 3 1.2 2.2 3.2 ` How exaclty to I regress my predictor variable x onto my data. I was thinking naively that I could just take the average of the data trials and get something like `x y 1 1.1 2 2.1 3 3.1 ` for which I could simply use `model <-lm(x ~ y)`. But I think I would introduce a bias by using the average. I was thinking instead to to create ordered pairs like (1,1.1), (1,2.0),...(3,2.2),(3,3.2). And use that for the regression but I'm not quite sure how to do this in a clean way in R. I was going to use data slicing and cbind a bunch of times, but I'm sure there is something better. **An example of what I'm getting at** The crux of my question is: how to do a regression model when you have data from mutiple runs of the same experiment. For example if I was looking at the temperature along a really thin metal wire (so its essentially 1-d), and then I mark off lenghts on this wire called x=0, x=1, x=2, ... x=total length. Then I take 10 temperature measuresments for each value of x. How do I create a regression model, which gives Temperarute as a function of length, i.e. T(x), given that I have 10 T values for each value of x?


title: Linear mixed effects model structure feedback

text: We are investigating the effects of a treatment (versus untreated control) on insect abundance. The study consists of 5 geographically separated study site replicates. Treatments were applied to 1 ~500 acre area at each location. Insects were collected from 14 randomly selected locations within the treatment area at each site, and were compared to samples collected from 14 randomly selected untreated reference area locations at each site. Samples were collected during 3 biweekly periods each year. So there are 840 observations for the entire study (2 years x 5 sites x 2 treatment levels/site x 14 sampling plots/treatment x 3 sampling periods/year). We are proposing to analyze these data with linear mixed effects models with treatment and years as the fixed effects, and site/sample location as random effects. Year is included as a fixed effect primarily because there are only 2 levels. I would appreciate feedback on whether either of the following models using lmer in the R lmer4 package are an appropriate starting point for analyzing these data:               m1 <- lmer(abundance ~ TRT*year + (period|site/sample))          Or               m2 <- lmer(abundance ~ TRT*year + (TRT|site) + (period|site/sample)      where `sample` is the 14 sampling locations within treatment level, and `period` is the 3 bi-weekly samples within each year.


title: Simultaneous confidence intervals for multinomial parameters, for small samples, many classes?

text: Suppose I have a sample from a multinomial distribution $\text{Multi}(p_1,\ldots,p_k)$, where $p_1+\dots+p_k=1$. I'd like simultaneous confidence intervals for the $p_i$'s. That is, for any $\alpha\in(0,1)$, I want to define $\ell_i$'s and $u_i$'s (random variables depending on the random sample) such that $$ P(p_1\in [\ell_1,u_1],\ldots,p_n\in[\ell_n,u_n])\ge 1-\alpha, $$ and obviously the closer the left hand side is to $1-\alpha$, the better. The special consideration is that the number of classes won't be too small (e.g. 10-100 classes), so I'm wondering if there's something significantly better than a Bonferroni approach. Thanks!


title: PCA's principle and available algorithms

text: If I understand correctly, PCA's principle is very simple:   1. Calculate data vectors' covariance matrix **C**.   2. Solve det( **C** - **e** *I) = 0, to find matrix **C** 's eigenvalues **e**.   3. Calculate matrix **C** 's eigenfunctions (from those eigenvalues). **FIRST:** Is this description correct? **SECOND:** Any algorithm for machine-solving of the polynomial equation det( **C** - **e** *I) = 0 ? I understand that this is a general math question (finding roots of a polynomial of degree **n** ). **THIRD:** Are there any **simple** implementations of PCA in C/C++?


title: Question about inverse-variance weighting

text: Suppose we want to make inference on an unobserved realization $x$ of a random variable $\tilde x$, which is normally distributed with mean $\mu_x$ and variance $\sigma^2_x$. Suppose there is another random variable $\tilde y$ (whose unobserved realization we'll similarly call $y$) that is normally distributed with mean $\mu_y$ and variance $\sigma^2_y$. Let $\sigma_{xy}$ be the covariance of $\tilde x$ and $\tilde y$. Now suppose we observe a signal on $x$, \begin{align}a=x+\tilde u,\end{align} where $\tilde u\sim\mathcal{N}(0,\phi_x^2)$, and a signal on $y$, \begin{align}b=y+\tilde v,\end{align} where $\tilde v\sim\mathcal{N}(0,\phi_y^2)$. Assume that $\tilde u$ and $\tilde v$ are independent. What is the distribution of $x$ conditional on $a$ and $b$? **What I know so far:** Using inverse-variance weighting, \begin{align}\mathbb{E}(x\,|\,a)=\frac{\frac{1}{\sigma_x^2}\mu_x+\frac{1}{\phi_x^2}a}{\frac{1}{\sigma_x^2}+\frac{1}{\phi_x^2}},\end{align} and \begin{align} \mathbb{V}\text{ar}(x\,|\,a)=\frac{1}{\frac{1}{\sigma_x^2}+\frac{1}{\phi_x^2}}. \end{align} Since $x$ and $y$ are jointly drawn, $b$ should carry some information about $x$. Other than realizing this, I'm stuck. Any help is appreciated!


title: Effect of switching response and explanatory variable in simple linear regression

text: Let's say that there exists some "true" relationship between $y$ and $x$ such that $y = ax + b + \epsilon$, where $a$ and $b$ are constants and $\epsilon$ is i.i.d normal noise. When I randomly generate data from that R code: `x <\- 1:100; y <\- ax + b + rnorm(length(x))` and then fit a model like `y ~ x`, I obviously get reasonably good estimates for $a$ and $b$. If I switch the role of the variables as in `(x ~ y)`, however, and then rewrite the result for $y$ to be a function of $x$, the resulting slope is always steeper (either more negative or more positive) than that estimated by the `y ~ x` regression. I'm trying to understand exactly why that is and would appreciate it if anyone could give me an intuition as to what's going on there.


title: Classification: model based and non model based

text: From discussion with Nick Cox, and some search, I was wondering what are the definitions of model-based classification and non-model-based classification? My understanding is that for model based classification, there is a statistical model for the joint distributions of the response and predictor. Some classification method is derived under some statistical model, such as linear discriminant analysis (LDA) assumes the conditional distribution of predictor given class is normal. Then what does non-model-based classification mean? Does it mean that the statistical model of the joint distributions of the response and predictor consists of any distributions? Can you given some classification methods for non-model-based classification? Thanks!


title: Rpart confidence interval

text: _Is there any way to compute 95% confidence intervals and prediction interval for the predicted value at each node?_ Thanks!!!


title: Comparing the distribution of several compositions given as percentages

text: I need to show that the proportion of different cell types in blood from two individuals are the same. Example: There are 6 cell types and                 Person 1 has the proportions: 0.07, 0.05, 0.4, 0.3, 0.15, 0.03   (Sums to 1)       Person 2 has the proportions: 0.07, 0.05, 0.3, 0.2, 0.35, 0.03   (Sums to 1)      Can anyone please suggest a test to show that the proportions in the two individuals are (or are not) significantly different?


title: When is it ok to merge data sets

text: There are two variables, $x$ and $y$, with $y$ depending on $x$ somewhat. For each occurrence of $x$, I record the value of $y$ so that I get a frequency distribution for $y$ given $x$. (The parentheses are to be read as (value, frequency) for $y$.)               x = 0: (19,2) (20,2) (35,5) (36,7)      x = 1: (14,1) (17,1)      x = 2: (2,1)      x = 3: (18,1) (19,2) (20,1)      x = 4: (27,1) (32,1)      x = 5: (19,1)      x = 7: (19,1)      x = 8: (19,1)      x = 9: (19,1)      x = 10: (19,1)      x = 11: (19,2)      x = 12: (19,3) (20,3) (25,1) (34,1)      x = 13: (19,5) (20,2)      x = 14: (19,3) (20,6) (21,1)      x = 15: (18,1) (19,316) (20,715) (21,4) (27,2) (33,1) (36,1) (37,1) (41,1)      x = 16: (18,2) (19,822) (20,462) (21,3) (22,1) (27,1) (30,1) (31,1) (34,1) (35,1) (39,1) (40,1) (42,1) (43,1)      x = 17: (18,1) (19,2) (20,1) (21,1)      x = 18: (17,1) (20,1)      x = 19: (16,1) (19,1) (20,1) (21,1)      x = 20: (15,1)      x = 23: (13,1)      x = 25: (10,1)      x = 26: (6,1)      x = 28: (6,2) (7,1)      x = 29: (6,1)      x = 32: (3,1)      x = 33: (2,1)      x = 37: (2,1)      As can be seen, $x$ is around 15 or 16 in the vast majority of cases, and the distribution for $y$ is about the same in those cases. So I would like to treat $x = 15$ and $x = 16$ as being the same. Furthermore, it looks as if I can lump more cases together; let's say $x = 12$ to $x = 17$ and conclude that if $12 \le x \le 17$ then we need no further information about $x$ to determine the distribution of $y$ (with some small error). What statistical procedures are relevant here? How can I numerically determine when distributions are similar enough? How can I aggregate the error?


title: How to reduce overfitting in linear regression

text: I am working with linear regression methods. The weakness of the method is the possibility of overfitting. So to reduce it, some papers use regularization. Are there other methods to reduce overfitting? Can we use a prior term to reduce overfitting? Given $D=\{(x_1,y_1);(x_2,y_2)...(x_n,y_n)\}$, the linear regression of the data $D$ is: $$H=wX+b$$ To reduce overfitting we add some regularization term. So the loss function is: $$J=\sum(h(x_i)-y_i)^2+\lambda_1\sum(w_i^2)$$ But finding $\lambda_1$ is so hard. Can we ignore it by using other terms to get more effective results? Thanks.


title: How do we define log-normal prior and a multivariate posterior log-likelihood in PyMC?

text: I am a newbie with pyMC and I am not still able to construct the structure of my MCMC with pyMC. I would like to establish a chain and I am confused how to define my parameters and log-likelihood function together. My chi-squared function is given by: ![enter image description here](http://i.stack.imgur.com/9o5rz.gif) where ![enter image description here](http://i.stack.imgur.com/wjSw9.gif) and ![enter image description here](http://i.stack.imgur.com/lXxgV.gif) are observational data and correspondence error respectively and ![enter image description here](http://i.stack.imgur.com/SaBGy.gif) is the model with four free parameter and the parameters are non-linear. The priors for `X` and `Y` are uniform but for `M` and `C` are given as following: ![enter image description here](http://i.stack.imgur.com/xy89G.gif) ; where the probability of `c` follows log-normal distribution while the expectation value of `c` is computed with the above formula and is the function of `M` and $\sigma$ is `0.09` if $M < 10^{15}$ otherwise $\sigma=0.06$: ![enter image description here](http://i.stack.imgur.com/OdM1n.gif) ![enter image description here](http://i.stack.imgur.com/Q779o.gif) for each `C` the parameter `z` is constant. I am wondering how I could define my likelihood for ![enter image description here](http://i.stack.imgur.com/XFaXo.gif) , and should it be referred as @Deterministic variable? Did I define `M` and `C` as priori information in a correct way or not? I will be grateful if somebody gives me some tips that how I can combine these parameters with given priors.               import pymc as pm     import numpy as np     import math     import random     from scipy.stats import expon          @pm.stochastic(dtype=np.float, observed=False, trace=True)     def Xpos(value=1900,x_l=1800,x_h=1950):         """The probable region of the position of halo centre"""         def logp(value,x_l,x_h):             if ((value>x_h) or (value<x_l)):            return -np.inf         else:            return -np.log(x_h-x_l+1)         def random(x_l,x_h):             return np.round((x_h-x_l)*random.random())+x_l          @pm.stochastic(dtype=np.float, observed=False, trace=True)     def Ypos(value=1750,y_l=1200,y_h=2000):         """The probable region of the position of halo centre"""         def logp(value,y_l,y_h):             if ((value>y_h) or (value<y_l)):            return -np.inf         else:            return -np.log(y_h-y_l+1)         def random(y_l,y_h):             return np.round((y_h-y_l)*random.random())+y_l               M=math.pow(10,15)*pm.Exponential('mass', beta=math.pow(10,15))          @pm.stochastic(dtype=np.float, observed=False, trace=True)     def concentration(value=4, zh, M200): #c parameter         """logp for concentration parameter"""         def logp(value=4.,zh, M):             if (value>0):            x = np.linspace(math.pow(10,13),math.pow(10,16),200 )            prob=expon.pdf(x,loc=0,scale=math.pow(10,15))            conc = [5.26/(1.+zh)*math.pow(x[i]/math.pow(10,14),-0.1) for i in range(len(x))]            mu_c=0            for i in range(len(x)):                mu_c+=prob[i]*conc[i]/sum(prob)            if (M < pow(10,15)):               tau=1./(0.09*0.09)            else:               tau=1./(0.06*0.06)                return  pm.lognormal_like(value, mu_c, tau)             else                return -np.inf         def random(mu_c,tau):             return np.random.lognormal(mu_c, tau, 1)


title: find correlated regions from time series

text: Given two time series, the problem is to find out correlated regions (by finding out the starting and ending points of those regions). For example, in the image: ![enter image description here](http://i.stack.imgur.com/c975u.png), the upper time series is a stock trend, the lower time series is market trend. In the region to the left of the red vertical bar, the two signals correlates quite well, meaning that the stock goes up and down together with the whole market. But in the region to the right of the red vertical bar, the stock trend has its own path, so in the right region, the two time series do not correlate well. The problem is how to find out such regions (there may be several regions in which the two time series correlates well).


title: Benefits of CART over ID3 algorithm

text: When building decision trees over a dataset that generates nodes with bad purity, is there any benefit of using the CART algorithm over the iterative dichotomizer 3 (ID3) algorithm?


title: Is it necessary to calculate Cronbach’s alpha for each sub-test?

text: I have a data set that includes two sub-tests, each measured using a Likert scale. I want to test the reliability. Is it necessary to calculate Cronbach’s Alpha for each sub-test? Or do I have to calculate just one Cronbach’s Alpha for all data?


title: I would like help calculating the probability of a simple problem

text: It's been a while since I've taken stat class so I would like help solving this problem. If the probability of a event occurring, like a person falling when they take a step, is 0.001%, then what is the probability of that person falling at least twice in a row in 1000 steps? Then what is the probability if 1 million people take 10 thousand steps that at least one person falls at least twice in a row?


title: testing differences between levels of a factor in a linear mixed model

text: I'm trying to wrap my head around using a linear mixed model and appropriate post-hoc tests to determine if there is a significant difference between various treatments in an experiment I inherited. (I'm really not sure if this question is more stats or more R, as it definitely contains elements of both.) The structure of the dataset is similar to that of the ChickWeight dataset in R, so I'll just use that as an example. (Four groups of chicks were put on different diets; each chick was weighed at multiple time points.) What I want to know is, does diet have a significant impact on the weight of the chicks? Specifically, is the effect of diet 2 significantly different from that of diet 3 or diet 4? (I don't much care about the effect of time, but need to account for it; the only way in which my data are actually different from the ChickWeight dataset is that all the chicks start off at approximately the same weight, and the starting numbers for my samples are significantly different in diet 2 vs 3 or 4.) So, I began my analysis by formulating a model (using the nlme package in R):               chicklme <- lme(weight~ Diet + Time + Diet:Time, random = ~1|Chick, data = ChickWeight, na.action = na.omit, method = "ML", correlation = corAR1(form = ~1 |Chick))      (With my dataset, the autoregressive correlation structure is definitely important.) To test whether Diet has a significant effect on weight, I compared my full model to the null model:               chicknull <- lme(weight~ Time + Diet:Time, random = ~1|Chick, data = ChickWeight, na.action = na.omit, method = "ML", correlation = corAR1(form = ~1 |Chick))     anova(chicknull, chicklme)      With the ChickWeights data, this difference actually isn't significant, but with my dataset, it is: lmeFull <- lme(weight ~ Diet + Time + Diet:Time, random = ~1|Chick, data = myData, na.action = na.omit, correlation = corAR1(form = ~1 |Chick), method = "ML") lmeNull <- lme(weight ~ Time + Diet:Time, random = ~1|Chick, data = myData, na.action = na.omit, correlation = corAR1(form = ~1 |Chick), method = "ML") > anova(lmeFull, lmeNull) Yields:                        Model df      AIC      BIC    logLik   Test  L.Ratio p-value     lmeFull     1 11 179.2363 218.9405 -78.61816                             lmeNull     2  8 205.1650 234.0408 -94.58249 1 vs 2 31.92865  <.0001      So, I conclude that Diet has a significant effect on weight. First, how do I report this? Bodo Winter has a very helpful tutorial (http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf), but I'm not entirely sure how to report this output (which differs from the anova output comparing two `lmer` models, but I can't use `lme4` because of the corAR1 business). Is it still doing a Chi squared test? What is the L.Ratio? `?anova` isn't very clear on this point. I know I can say "diet had a significant effect on weight," but aside from the p-value, I'm not sure what statistical info I need to put into the description of the result (when I'm writing it up for publication). Next, I have the issue of comparing the different diets. I've been using `glht` from the `multcomp` package, but I'm still trying to make sure I understand what I am asking it to test. If I then do a Tukey HSD test on diet 2 vs 3 and diet 2 vs 4 like so               lmeFull <- update(lmeFull, method = "REML")     summary(glht(lmeFull, lnfct = mcp(Diet = c("2-3=0", "2-4=0"))))      am I correct in thinking that this tests the effect of diet on weight only on the first level of Time? Since the chicks all started out at more or less the same weight, this is not very helpful; I know I can center the time variable at any time I want, and for this analysis, the last measurement makes most sense, but I have no idea how to do this in R, since it appears to take the lowest value of the numerical factor Time as the reference level. And, if the post-hoc Tukey test does only test differences at a particular time point, and not differences overall, is there an appropriate way to try to do that? Or do I just pick a time point, and analyze that? I know that the Time:Diet interaction (significant in my dataset) tells me that the change in weight over time differs between diets, but if what I want to know is "do the chicks on diet 2 weigh less than the chicks on diets 3 and 4?", what is the most appropriate test to use? Can I just remove the Time factor from the model and test the weights from all time points? This seems rather silly, as it's throwing so much of the data away. I've been reading and reading to try to understand, and I think there is something about the stats that I am fundamentally missing. If anyone can explain this to me in plain English, I would truly appreciate it!


title: Difference between training, test and holdout set data mining model building

text: What is the difference between training, test, and holdout sets? I know these concepts, just want to ensure that I have understood correctly. Training set is something that we have as of now. We will remove subset from it and removed subset will be called holdout set. We will build models using remaining data (what remains after removing holdout set) and the holdout set is used to finalized estimates of tuning parameters **_(step 1)_** Then we will build a final model on the entire Training set (including holdout set). Tuning parameters values are same as that we got from step 1. Test data is something that we get in future. We don't know their Y/dependent variable value and we predict it using our model.


title: Feature selection with Random Forests

text: I have a dataset with mostly financial variables (120 features, 4k examples) which are mostly highly correlated and very noisy (technical indicators, for example) so I would like to select about max 20-30 for later use with model training (binary classification - increase / decrease). I was thinking about using random forests for feature ranking. Is it a good idea to use them recursively? For example, let's say in the first round I drop the worst 20%, second too and so on until I get the desired number of features. Should I use cross-validation with RF? (It's intuitive for me not to use CV because that's pretty much what RF does already.) Also if I go with random forests should I use them as classifiers for the binary or regressor for the actual increase / decrease to get feature importances? By the way, the models I would like to try after feature selection are: SVM, neural nets, locally weighted regressions, and random forest. I'm mainly working in Python.


title: Thoughts on best method for testing for differences in slopes of lines (time series data)?

text: I have 4 sets of 50 year data (one value per year for each of the 4 different categories)- and would like to test for differences in the slopes of the lines. I think ANCOVA is appropriate- though maybe not as these are time series data? And ...how to do post hoc pairwise comparisons? The data are estimates of parasite production for each year under 4 different climate scenarios ranging from hot and dry to cold and wet (parasite production is directly dependent on temperature).


title: Unit root test confusion

text: For my time-series regression, I am regressing the difference in variable x on a difference in variable y. Before proceeding, I want to check for stationarity of my variables. Regressing d.x (difference in x) over L.d.x (lag of the difference in x), I get a coefficient of $-0.9995261$ and a 95% interval: $[-1.000654, -0.9983979]$. As far as I remember, the abs value of the coefficient (and the values in the interval) must be less than 1 for us to be sure of the stationarity, which is not the case here. So I conducted a Dickey-Fuller test (all of this in Stata) and got the following results:               Interpolated Dickey-Fuller              Test Statistic  |  1% Critical Value |  5% Critical Value | 10% Critical Value           Z(t)       -3474.557   |         -3.430     |     -2.860         |    -2.570          MacKinnon approximate p-value for Z(t) = 0.0000      I am confused because they very strongly suggest that there is no unit root, whereas I expected a close call. Is this because I have ~5000 observations? Or am I wrong, and the condition for the unit root is a coefficient in excess of 1, not its absolute value in excess of one?


title: Time series classification & hierarchal clustering

text: I have a need to build a classifier for time series data. I've made progress but would like advise on my current challenge (see below). My research led me to KNN and Dynamic Time Warping as the distance measure. This has performed very well on a small dataset. Unfortunately, it is computationally intensive. Paper: http://machinelearning.wustl.edu/mlpapers/paper_files/icml2006_XiKSWR06.pdf So, I want to select observations to be included in my training set for KNN. It is important that these observations are representative of the overall data. I would like feedback and thoughts on the following idea:   * Generate a pairwise distance matrix using DTW    * Perform some for of clustering (perhaps hierarchal) to identify time series that are similar to one another   * Randomly sample observations from these clusters to generate my training set Questions:   * Does this approach make sense?   * Can you advise on a clustering method?   * Any red flag? Advice? Much appreciated.


title: Enormous coefficients in logistic regression - what does it mean and what to do?

text: I get enormous coefficients during logistic regression, see coefficients with `krajULKV`:               > summary(m5)          Call:     glm(formula = cbind(ml, ad) ~ rok + obdobi + kraj + resid_usili2 +          rok:obdobi + rok:kraj + obdobi:kraj + kraj:resid_usili2 +          rok:obdobi:kraj, family = "quasibinomial")          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.7796  -1.0958  -0.3101   1.0034   2.8370            Coefficients:                                   Estimate     Std. Error t value Pr(>|t|)        (Intercept)                 -486.72087      664.71911  -0.732  0.46424        rok                            0.24232        0.33114   0.732  0.46452        obdobinehn                  3400.43703     1354.14874   2.511  0.01223 *      krajJHC                      786.22409      708.50291   1.110  0.26746        krajJHM                      511.85538      823.03038   0.622  0.53417        krajLBK                      -23.94180     2388.86316  -0.010  0.99201        krajMSK                     1281.88767      955.09736   1.342  0.17992        krajOLK                     -175.19425     1255.82946  -0.140  0.88909        krajPAK                      349.76438     1071.03364   0.327  0.74408        krajPLK                    -1335.73206     1534.09899  -0.871  0.38418        krajSTC                      868.99157      692.30426   1.255  0.20976        krajULKV                  245661.86828 17496742.31677   0.014  0.98880        krajVYS                     3341.76686     1314.77140   2.542  0.01121 *      krajZLK                     3950.75617     2922.25220   1.352  0.17676        resid_usili2                  -1.44719        0.89315  -1.620  0.10555        rok:obdobinehn                -1.69479        0.67462  -2.512  0.01219 *      rok:krajJHC                   -0.39108        0.35295  -1.108  0.26817        rok:krajJHM                   -0.25481        0.40997  -0.622  0.53443        rok:krajLBK                    0.01621        1.19155   0.014  0.98915        rok:krajMSK                   -0.63985        0.47592  -1.344  0.17917        rok:krajOLK                    0.08714        0.62545   0.139  0.88923        rok:krajPAK                   -0.17419        0.53344  -0.327  0.74410        rok:krajPLK                    0.66539        0.76383   0.871  0.38394        rok:krajSTC                   -0.43292        0.34490  -1.255  0.20976        rok:krajULKV                -122.01076     8704.03367  -0.014  0.98882        rok:krajVYS                   -1.66391        0.65468  -2.542  0.01122 *      rok:krajZLK                   -1.96718        1.45474  -1.352  0.17667        obdobinehn:krajJHC         -3623.86807     1385.86009  -2.615  0.00909 **     obdobinehn:krajJHM         -3220.08906     1458.83842  -2.207  0.02757 *      obdobinehn:krajLBK         -1051.07131     3434.11845  -0.306  0.75963        obdobinehn:krajMSK         -6415.65781     1978.30260  -3.243  0.00123 **     obdobinehn:krajOLK         -2427.66591     1777.51914  -1.366  0.17239        obdobinehn:krajPAK         -3111.45312     1623.59145  -1.916  0.05566 .      obdobinehn:krajPLK         -1800.26258     2065.74461  -0.871  0.38375        obdobinehn:krajSTC         -4409.45624     1379.64196  -3.196  0.00145 **     obdobinehn:krajULKV      -187832.68360 16454272.74951  -0.011  0.99089        obdobinehn:krajVYS         -5445.51446     1791.38012  -3.040  0.00244 **     obdobinehn:krajZLK         -6216.43343     3167.49836  -1.963  0.05003 .      krajJHC:resid_usili2           1.60474        0.98554   1.628  0.10385        krajJHM:resid_usili2           1.57822        1.04518   1.510  0.13143        krajLBK:resid_usili2          11.53462       13.40012   0.861  0.38961        krajMSK:resid_usili2          -1.33600        1.55241  -0.861  0.38971        krajOLK:resid_usili2           0.07296        1.27034   0.057  0.95421        krajPAK:resid_usili2           1.35880        1.23033   1.104  0.26974        krajPLK:resid_usili2           1.90189        1.41163   1.347  0.17826        krajSTC:resid_usili2           2.05237        0.95972   2.139  0.03277 *      krajULKV:resid_usili2        599.79215    20568.86123   0.029  0.97674        krajVYS:resid_usili2           3.03834        1.16464   2.609  0.00925 **     krajZLK:resid_usili2           1.18574        1.11024   1.068  0.28583        rok:obdobinehn:krajJHC         1.80611        0.69042   2.616  0.00906 **     rok:obdobinehn:krajJHM         1.60475        0.72676   2.208  0.02751 *      rok:obdobinehn:krajLBK         0.52268        1.71244   0.305  0.76027        rok:obdobinehn:krajMSK         3.19712        0.98564   3.244  0.00123 **     rok:obdobinehn:krajOLK         1.21012        0.88541   1.367  0.17208        rok:obdobinehn:krajPAK         1.55034        0.80886   1.917  0.05563 .      rok:obdobinehn:krajPLK         0.89718        1.02893   0.872  0.38349        rok:obdobinehn:krajSTC         2.19742        0.68732   3.197  0.00144 **     rok:obdobinehn:krajULKV       93.43130     8189.24994   0.011  0.99090        rok:obdobinehn:krajVYS         2.71357        0.89236   3.041  0.00243 **     rok:obdobinehn:krajZLK         3.09624        1.57711   1.963  0.04996 *      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for quasibinomial family taken to be 1.258421)              Null deviance: 1518.0  on 878  degrees of freedom     Residual deviance: 1228.6  on 819  degrees of freedom       (465 observations deleted due to missingness)     AIC: NA          Number of Fisher Scoring iterations: 18      **What does this mean??** Does it mean some multicollinearity, like @Scortchi meantioned in this discussion? Or does this mean overfitting? How to detect the problem? **What shall I do now?** I tried to remove some variables. This helps a bit but not so much:               > m6 <- update(m5, ~.- kraj:resid_usili2)     > m7 <- update(m6, ~.- resid_usili2)     > summary(m7)          Call:     glm(formula = cbind(ml, ad) ~ rok + obdobi + kraj + rok:obdobi +          rok:kraj + obdobi:kraj + rok:obdobi:kraj, family = "quasibinomial")          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.9098  -1.1931  -0.2274   1.0529   3.1283            Coefficients:                                Estimate  Std. Error t value Pr(>|t|)     (Intercept)              -118.95199   476.34698  -0.250    0.803     rok                         0.05971     0.23718   0.252    0.801     obdobinehn                412.69412   646.95083   0.638    0.524     krajJHC                   447.69791   498.45358   0.898    0.369     krajJHM                   -62.92516   525.85737  -0.120    0.905     krajLBK                   677.73239  1595.20024   0.425    0.671     krajMSK                   278.24639   621.32312   0.448    0.654     krajOLK                  -705.97832   782.53474  -0.902    0.367     krajPAK                   387.96543   608.98961   0.637    0.524     krajPLK                  -653.68419   782.20737  -0.836    0.403     krajSTC                  -114.34822   489.06318  -0.234    0.815     krajULKV                -2117.64674  1797.75836  -1.178    0.239     krajVYS                   884.74411   681.05324   1.299    0.194     krajZLK                  -997.77613   925.93280  -1.078    0.281     rok:obdobinehn             -0.20602     0.32211  -0.640    0.523     rok:krajJHC                -0.22303     0.24819  -0.899    0.369     rok:krajJHM                 0.03092     0.26180   0.118    0.906     rok:krajLBK                -0.33909     0.79438  -0.427    0.670     rok:krajMSK                -0.13889     0.30935  -0.449    0.654     rok:krajOLK                 0.35102     0.38943   0.901    0.368     rok:krajPAK                -0.19392     0.30323  -0.640    0.523     rok:krajPLK                 0.32463     0.38937   0.834    0.405     rok:krajSTC                 0.05677     0.24351   0.233    0.816     rok:krajULKV                1.05287     0.89453   1.177    0.239     rok:krajVYS                -0.44149     0.33911  -1.302    0.193     rok:krajZLK                 0.49612     0.46081   1.077    0.282     obdobinehn:krajJHC       -776.31258   672.68911  -1.154    0.249     obdobinehn:krajJHM       -267.78650   700.38741  -0.382    0.702     obdobinehn:krajLBK      -1246.67321  1760.37329  -0.708    0.479     obdobinehn:krajMSK       -383.77613   858.81391  -0.447    0.655     obdobinehn:krajOLK        -96.72334   947.75189  -0.102    0.919     obdobinehn:krajPAK       -540.25140   827.13134  -0.653    0.514     obdobinehn:krajPLK       -517.49161  1124.63474  -0.460    0.645     obdobinehn:krajSTC       -683.81160   672.66674  -1.017    0.310     obdobinehn:krajULKV      2344.32314  2073.98366   1.130    0.259     obdobinehn:krajVYS       -795.62043   917.80551  -0.867    0.386     obdobinehn:krajZLK        618.33075  1093.37768   0.566    0.572     rok:obdobinehn:krajJHC      0.38725     0.33493   1.156    0.248     rok:obdobinehn:krajJHM      0.13374     0.34870   0.384    0.701     rok:obdobinehn:krajLBK      0.62237     0.87662   0.710    0.478     rok:obdobinehn:krajMSK      0.19114     0.42758   0.447    0.655     rok:obdobinehn:krajOLK      0.04842     0.47171   0.103    0.918     rok:obdobinehn:krajPAK      0.26922     0.41184   0.654    0.513     rok:obdobinehn:krajPLK      0.25790     0.55986   0.461    0.645     rok:obdobinehn:krajSTC      0.34078     0.33492   1.017    0.309     rok:obdobinehn:krajULKV    -1.16571     1.03236  -1.129    0.259     rok:obdobinehn:krajVYS      0.39675     0.45704   0.868    0.386     rok:obdobinehn:krajZLK     -0.30732     0.54422  -0.565    0.572          (Dispersion parameter for quasibinomial family taken to be 1.313286)              Null deviance: 2396.8  on 1343  degrees of freedom     Residual deviance: 2110.3  on 1296  degrees of freedom     AIC: NA          Number of Fisher Scoring iterations: 5      **EDIT:** As proposed by Scortchi, I tried to use VIF and I also get enormous values. What does this mean? See:               > require(HH)     > vif(cbind(ml, ad) ~ rok + obdobi + kraj + resid_usili2 +      +         rok:obdobi + rok:kraj + obdobi:kraj + kraj:resid_usili2 +      +         rok:obdobi:kraj)                         rok              obdobinehn                 krajJHC                 krajJHM                    50.281603         45075363.969712         15194580.406796         11362184.620230                      krajLBK                 krajMSK                 krajOLK                 krajPAK               7567915.376763          5228018.864051         17105623.986998         10944471.683601     [... cut out ...]


title: Means in R: How to find this when the exact list is not known

text: This is a simple question, not really related to coding but more so to some math that I am missing in my brain. I know: An object A has a length of 5 and a mean of 4. An object B has a length of 10 and a mean of 7. I make a new object, C, where C <- c(A, 1, B). What is the mean of C?


title: Weighted Linear Regression R

text: Can anyone expalin to me in simple terms what happens when we use weights in `regsubsets` or `lm` in R? What effect do weights have on a linear regression? for example :               Model1<-lm(Ozone~Solar.R,data=airquality)     summary(Model1)     #Coefficients:     #            Estimate Std. Error t value Pr(>|t|)         #(Intercept) 18.59873    6.74790   2.756 0.006856 **      #Solar.R      0.12717    0.03278   3.880 0.000179 ***     Model1<-lm(Ozone~Solar.R,data=airquality,weights=(2*seq(nrow(airquality),1,-1)))     summary(Model1)     #Coefficients:     #            Estimate Std. Error t value Pr(>|t|)         #(Intercept) 18.57106    6.26067   2.966 0.003704 **      #Solar.R      0.10824    0.02927   3.699 0.000341 ***      please explain the changes in intercepts and slope.


title: Different possible parameter combinations to obtain a standardized generalized hyperbolic distribution?

text: Consider the generalized hyperbolic distribution given by (from wikipedia): ![gh](http://i.stack.imgur.com/GYWx0.png) So I now wanted to derive the standardized version, so mean zero variance one. I wanted to do the following:   1. Set the mean to zero and the variance equal to one.   2. Rewrite the mean equation, such that  \begin{align*}\beta=-\mu\frac{\gamma K_\lambda (\delta \gamma)}{\delta K_{\lambda+1}(\delta\gamma)}\end{align*}   3. Insert this expression into the variance equation.   4. Rewrite the equation, such that $\mu =.....$   5. Insert this expression into the pdf. This gives the formula for the standardized ghyp, which now only has the paramger $\lambda$, $\alpha$, $\beta$, $\delta$, since $\mu$ is fixed. My question is now: Can I do this, I mean, is this correct? My second important question is: I think, there are different ways to express the point, that the mean is zero and the variance equal to one. So instead of solving for $\beta$ and then $\mu$, I could have solved for different combinations of the parameters to get the mean zero and the variance equal to one. So my solution is just one possibility. Consider the following simplified example (arbitrary equations, just for illustration): \begin{align*} \text{GH}= \mu + \alpha + 3 \beta + 4 \delta \end{align*} The mean e.g. is given by and set to zero: \begin{align*} \text{mean}= 0 = \mu + \delta \beta \end{align*} The variance e.g. is given by and set to one: \begin{align*} \text{variance} = 1 = \beta \delta +\delta /(\alpha-\beta) \end{align*} There are now different ways of rewriting and solving. I can compute a solution, in which the final equation depends on $\alpha$, $\beta$, and $\delta$ and $\mu$ fixed, or, there is e.g. a solution in which the final equation depends only on $\mu$, and $\delta$ and $\alpha$ and $\delta$ is fixed. So both would fulfill the requirement of mean zero and variance one. I think by going the way I proposed in case of the GHYP above, there is the same problem, that this is only one solution, but there exist other combinations, in which the mean zero and variance one equation is still fulfilled. So: What would you do in such a case?


title: Confusion with Augmented Dickey Fuller test

text: I am working on the data set `electricity` available in R package `tsa`. My aim is to find out if an `arima` model will be appropriate for this data and eventually fit it. So I proceeded as follows:      1st: Plot the time series which resulted if the following graph: ![ts plot1](http://i.stack.imgur.com/2t3Qk.png)      2nd: I wanted to take log of `electricity` to stabilize variance and afterward differenced the series as appropriate, but just before doing so, I tested for stationarity on the original data set using the `adf` (Augmented Dickey Fuller) test and surprisingly, it resulted as follows:      ### Code and Results:               adf.test(electricity)                       Augmented Dickey-Fuller Test     data:  electricity      Dickey-Fuller = -9.6336, Lag order = 7, p-value = 0.01      alternative hypothesis: stationary     Warning message: In adf.test(electricity) : p-value smaller than printed p-value      Well, as per my beginner's notion of time series, I suppose it means that the data is stationary (small p-value, reject null hypothesis of non- stationarity). But looking at the ts plot, I find no way that this can be stationary. Does anyone has a valid explanation for this?


title: Confusion with Chi-Square test

text: I have an extremely simple question regarding chi-square test. So far I have found two different formulas for it: $$ \chi^2=\sum\frac{(x_o-x_e)^2}{x_e} $$ and $$ \chi^2=\sum\frac{(x_o-x_e)^2}{\sigma^2} $$ $x_o$ is the observed data while $x_e$ is the expected data. Clearly the two formulas do not return the same result because $x_e$ != $\sigma^2$. Can someone explain to me why there are two different formulas for the chi-square test and when I should use each one of them? Thank you.


title: Positive correlations to dependent variable, but negative coefficients

text: First of all, sorry for the huge pictures, but I'm in desperate need of some input and help on the results following a study. I'm trying to interpret the results from my study, but I can't quite figure out what causes some of the positive correlations to turn into negative coefficients in the regression analysis. ![Correlations](http://i60.tinypic.com/snnsb4.jpg) As you can see there are positive correlations between "atmosfære", "kampkvalitet", "fotball" and the dependent variable "intensjon". What I don't understand is how these three independent variables end up with negative coefficients in the following regression analysis. Considering VIF-values are low, which indicates low multicollinearity, I can't seem to wrap my head around this. ![Regression](http://i62.tinypic.com/ibi652.jpg) Some thoughts and input that might help me along the way? Much obliged!


title: SUR estimation and Heckman selection model with panel data on Stata?

text: I'm working with unbalanced panel data using time and firms as IDs and would like to find out how to test for correlation between two panel equations that may be seemingly unrelated. The second question is about any package that allows the use of a Heckman selection model for panel data in Stata.


title: Smoothing Spline Example

text: I am learning the smoothing spline method. I saw that smoothing spline is a penalty term to reduce overfitting in linear regression. Given dataset {$(x_1,y_1),(x_2,y_2)..(x_n,y_n)$}So the formular such as: $$RSS=\sum(y_i-f(x_i))^2+\lambda\int((f(t)'')^2dt$$ Assume it is linear case so $$f(x_i)=ax_i+b$$ $$f(x_i)''=0$$ Is it correct. Could you explaint help me how to find second term in RSS ($\lambda\int((f(t)'')^2dt$) in linear regression case? Or give me one example?Thank you so much


title: Logistic regression: Strange standard errors from glm() in R

text: To my surprise I found that standard errors and thus Wald confidence intervals became smaller when I removed the intercept from a simple logistic regression model, using `glm()` and R.               # load an object named "my.df" to the global enviroment     load(url("http://hansekbrand.se/code/test.RData"))     # fit a model with intercept to data     my.fit <- glm(deprived.of.education ~ religion, data = my.df, family = binomial("logit"))     # fit a model without any intercept to data     my.fit.without.intercept <- glm(deprived.of.education ~ 0 + religion, data = my.df, family = binomial("logit"))          # inspect the first fit     summary(my.fit)$coefficients     #                        Estimate Std. Error   z value     Pr(>|z|)     # (Intercept)          -2.8718056 0.03175130 -90.44687 0.000000e+00     # religionChristianity  0.4934891 0.03234887  15.25522 1.519805e-52     # religionHinduism      0.5257316 0.03376535  15.57015 1.161317e-54     # religionIslam         1.5734832 0.03231692  48.68914 0.000000e+00     # religionNonreligious  1.5975456 0.03555164  44.93592 0.000000e+00          # inspect the second fit     summary(my.fit.without.intercept)$coefficients     #                       Estimate  Std. Error    z value Pr(>|z|)     # religionBuddhism     -2.871806 0.031751299  -90.44687        0     # religionChristianity -2.378317 0.006189045 -384.27842        0     # religionHinduism     -2.346074 0.011487113 -204.23530        0     # religionIslam        -1.298322 0.006019850 -215.67354        0     # religionNonreligious -1.274260 0.015992939  -79.67642        0      I understand why the _z_ values are different, because the null hypotheses in the two cases are different. In the first case, with the intercept, the null is "same as the reference category", while without the intercept, the null becomes "zero".   But I do not understand the large difference in standard errors between the two models. Without the intercept, the standard errors seem to vary with _n_ of each level, i.e. there are many cases of "Christianity" and "Islam", and they have small standard errors, but with the intercept, there is essentially no variation in the standard errors. Could someone please explain the reason for the differences in the magnitude of the standard errors between the two models? I would like to calculate probabilities and confidence intervals around them, and I have done so using the estimates from the first model. If I would do that with the estimates from the second model, the confidence intervals would be much smaller, but would they be reliable?


title: Measurment error for two variables

text: I am interested in estimating the effect of security S on crime C in a given city over time (eight quarters) for twenty cities, so it's panel data. The problem is, instead of actual security spending I have two imperfect measures of S, call them S1 and S2. These were collected independently from each other but in both measures some random mistakes caused mis-measurment. I know that if I use either of the two in order to proxy S, my regression coefficient will be biased towards zero because of attenuation. Now I have two questions: why does attenuation bias lead to downward biased coefficients? Can I use S1 and S2 somehow to estimate the causal effect of S on C? Thanks


title: Find the midpoint of S-curves with R

text: I have a bunch of curves that look like S-curves _(stored in a db with Point format)_ and I am interested in finding this midpoint between both horizontal asymptotes for each of them. So, if **X-axis** represents **time** and **Y-axis** number of **occurrences of a phenomenon** , how can I project this midpoint over X-axis to know the date it happened? I would like to do this operation in R software, please, can you help me with some guidelines about how to do it? Thanks!


title: Getting the units right for the Pareto distribution of wealth: x = people, dollars, dollars per person?

text: When people talk about the 80-20 rule in the context of wealth, it is usually expressed, verbally, by stating that the 20 percent of the people with the highest wealth get 80 percent of the wealth, and the top 20 percent of the richest 20 get 80 percent of that, or 64% of total wealth, and so forth. Then there is hand-waving about how the Pareto distribution represents this constant proportionality. However, when you graph a Pareto distribution in income, the x value, income, has to be measured either in dollars or in dollars per person. This implies that the pdf is measured in 1/dollars or in persons per dollar, respectively, so that the area of the integral can be a unitless probability. If we measure income in dollars, the interpretation is straightforward. The CDF F(x) gives the probability that an arbitrary dollar will be found in the pocket of someone with an income below x. But this relates dollars to dollars, or a unitless proportion to dollars of income per person. Suppose I want to know about the income of groups of people, as in the popular explanation with which I opened. I have not been able to make either the dollars to (1/dollars) units or the dollars per person to persons per dollar units come out in a way that gives me a percent of people, or of people in a particular income range, instead of a percent of dollars in all dollars. Is there someone who could give a clear explanation of how you set the units up so as to get the answer in terms of a ratio between the percent of people and the percent of wealth? Assuming that the value of alpha is set at the level that give you the 80/20 rule exactly, I'd like to see what, if this were chemistry, I'd call the stoichiometry: an arrangement of the problem so that the units on the left hand side are equal to the units on the RHS, and those units come out relating a proportion of people to a proportion of wealth, e.g. showing that the top 1/25 of the people get 64% of the wealth.


title: Can we conclude structural break from the following results?

text: The following figure shows the log-likelihood of individual observations for different subsets of a data set - the horizontal line shows the average log- likelihood: ![figure one](http://i.stack.imgur.com/14O5r.png) So, the first subset includes all the data points observed up to t=7, the second all data points up to t=33, and so forth. For each subset, I estimate a parameter $\theta$, and then I calculate the likelihood of each observation given this parameter estimate. This is what the figure shows. It can be seen that on average the likelihood of the observed data (shown as a horizontal line) tends to decrease as the number of observations grows. This pattern can be confirmed in the following plot, which shows how the average likelihood of observations changes as the number of observations included grows, for different data sets: ![enter image description here](http://i.stack.imgur.com/w3c8V.png) In the light of this, would you agree that these results suggest that the $\theta$ parameter is not equal for all $t$? If it was equal, I would expect the same average likelihood regardless of the number of observations included.


title: Least Squares Regression in Stata

text: Can anyone explain how Stata computes the least squares parameters with a single explanatory variable - what algorithm Stats uses to compute the parameters. Actually I want to understand the computational algorithm of negative binomial regression model.


title: Using panel OLS coefficients for forecasting - gretl

text: I conducted a panel OLS regression using gretl, on a balanced panel with around 25 zones and 10 time periods. I would like to use the model to forecast future values of y (dependant variable), given known future changes of x (explanatory variable). I would like to ask: how to deal with (a) the constant, (b) the unit dummies; (c) the time dummies when it comes to forecasting. Do i have to include any/all of the above in computing y-hat, or simply use the coefficient for the x term?


title: is psych::principal function still PCA when using rotation?

text: I have tried to reproduce some research (using PCA) from SPSS in R. In my experience, `principal()` function from package `psych` was the only function that came close (or if my memory serves me right, dead on) to match the output. To match the same results as in SPSS, I had to use parameter `principal(..., rotate = "varimax")`. I have seen papers talk about how they did PCA, but based on the output of SPSS and use of rotation, it sounds more like Factor analysis. Question: Is PCA, even after rotation (using `varimax`), still PCA? I was under the impression that this might be in fact Factor analysis... In case it's not, what details am I missing?


title: Open Source statistical textbooks?

text: There have been a few questions about statistical textbooks, such as the question Free statistical textbooks. However, I am looking for textbooks that are Open Source, for example, having an Creative Commons license. The reason is that in course material in other domains, you still want to include some text about basic statistics. In this case, it would be interesting to reuse existing material, instead of rewriting that material. Therefore, what Open Source textbooks on statistics (and perhaps machine learning) are available?


title: JAGS - problem with taking submatrices

text: I'm trying to use JAGS for clustering mixtures of multivariate normal distributions. In order to model different covariance structures in each cluster I wanted to use one big (C*D)xD matrix (C - number of clusters, D - number of dimensions) which would represent all of the covariance matrices stacked one upon another. The problem is, that while JAGS is perfectly happy with sampling such submatrices from Wishart distribution, it will not allow me to subsequently use them as parameters for the normal distribution (or at least so I suspect). Here is my (not so) minimal working example:               library(rjags)     library(MASS)          clusters.bug <-     "model {         # N - observations         # D - dimensions         # C - clusters              for (i in 1:N) {             # JAGS has some problem with the following indexing here:             x[i,1:D] ~ dmnorm(mu[c[i],1:D], Omega[(D*(c[i]-1)+1):(D*c[i]),1:D])             # x[i,1:D] ~ dmnorm(mu[c[i],1:D], Omega0)             c[i] ~ dcat(p)         }              for (j in 1:C) {             mu[j,1:D] ~ dmnorm(mu0, Omega0)             Omega[(D*(j-1)+1):(D*j),1:D] ~ dwish(Omega0, D)         }              for (k in 1:D) {             mu0[k] <- 0             for (l in 1:D) {                 Omega0[k,l] <- ifelse(k == l, 1, 0)             }         }              p ~ ddirch(alpha)     }"          ## Data     N <- 1000     D <- 2     C <- 3          p <- c(0.2, 0.3, 0.5)     mu <- matrix(c(1,2, 3,4, 5,6), ncol = 2, byrow = TRUE)     cv <- matrix(c(1,0, 0,1, 1,.5, .5,2, 2,.3, .3,1), ncol = 2, byrow = TRUE) / 10          x <- NULL     for(i in 1:N) {         k <- sample(1:C, 1, prob = p)         x <- rbind(x, mvrnorm(1, mu[k,], cv[(D*(k-1)+1):(D*k),]))     }     contour(kde2d(x[,1], x[,2]))          alpha <- rep(1, C)     Omega <- NULL     for(i in 1:C) {         Omega <- rbind(Omega, diag(D))     }          ## JAGS     jags <- jags.model(textConnection(clusters.bug),                        data = list('x' = x,                                    'N' = N,                                    'D' = D,                                    'C' = C,                                    'alpha' = alpha),                        n.chains = 1,                        n.adapt = 100)      If you comment out:               x[i,1:D] ~ dmnorm(mu[c[i],1:D], Omega[(D*(c[i]-1)+1):(D*c[i]),1:D])      and uncomment the next line, everything is ok.


title: Independent and serial correlation

text: I have a set of 76 annual occurrences of $x$ and $y$. When I run a linear regression, $y$ is significantly correlated to $x$. The Durbin-Watson test shows $y$ is also serially correlated. My goal is to predict $y$ from a certain future $x$ and to determine the prediction interval. I believe the proper point estimate equation is: $$Y_t = X_tB + pY_{t-1} – pX_{t-1}B + V_t.$$ Is this correct? $X_t$, $Y_{t-1}$, and $X_{t-1}$ are known. What are $B$, $p$, and $V_t$ and how can I calculate them for a particular $Y_t$ prediction? How do I calculate the prediction interval for $Y_t$? If not too complex, how does equation change if serial correlation is second order?


title: Implications of lower-bounded total variation distance on hypothesis testing

text: Let $\{X_i\}_n$ be a sequence of $n$ random variables independently and identically drawn from either $P$ or $Q$. Thus the sequence $\{X_i\}_n$ has a product distribution, which is either $P^n$ or $Q^n$. Assume that $P$ and $Q_\theta$ have the same support $\Omega$. Suppose that the total variation distance between $P^n$ and $Q^n$ is lower- bounded: $$\delta(P^n,Q_\theta^n)=\frac{1}{2}\int_{\Omega}|f_P^n(x)-f_{Q}^n(x)|dx\geq 1-\epsilon_n$$ (here $f_P^n(\cdot)$ and $f_{Q}^n(\cdot)$ are respective densities of $P^n$ and $Q^n$) Now suppose that $\epsilon_n$ converges to zero as $n$ gets large, i.e. $\lim_{n\rightarrow \infty}\epsilon_n=0$. We can also assume uniform convergence, if that helps. Basically, by choosing $n$ large enough, we can obtain $\delta(P^n,Q^n)$ arbitrarily close to unity. Does this imply that there **exists** a sequence of hypothesis tests $\phi_n(\{X_i\}_n)$ that classifies $\{X_i\}_n$ as drawn either from $P^n$ or $Q^n$ with the sum of error probabilities $\alpha_n+\beta_n\leq\epsilon_n$, where $\alpha_n$ is the probability of mis-classifying $\{X_i\}_n$ that was drawn from $P^n$ as drawn from $Q^n$, and $\beta_n$ is the probability of mis- classifying $\{X_i\}_n$ that was drawn from $Q^n$ as drawn from $P^n$? If so, is there a well-defined procedure on determining how to construct such a sequence of tests? I know that $1-\delta(P^n,Q^n)$ provides a lower bound on the sum $\alpha_n+\beta_n$ for any hypothesis test, but I am not sure what happens if $\delta(P^n,Q^n)$ is lower-bounded close to unity.


title: Variance-covariance matrix

text: $\DeclareMathOperator{\var}{Var}$ How to compute prediction bands for non- linear regression? In the above link, you have mentioned about the variance-covariance matrix of the estimates. What is the explicit formula for computing the terms of this variance-covariance matrix? If suppose, $f = f(x;p1,p2,p3)$ with $x$ as the independent variable and $P=(p1,p2,p3)$ are the nonlinear parameters of a function which fits the data in a nonlinear regression problem, what is $\var(P)$? Is the Hessian equal to the one presented in this link? https://www8.cs.umu.se/kurser/5DA001/HT07/lectures/lsq-handouts.pdf or http://www.math.vt.edu/people/dlr/m2k_svb11_hesian.pdf Or, are they the same? Can you explain the difference? Do the inverse of the above Hessians mentioned in the above links estimate the variance-covariance matrix? Which one of the Hessian formula must be taken into consideration for evaluating the Hessian you meant in the first link mentioned above? Rizwan


title: Confidence interval for the difference between $X_1$ and $X_2$ or $X_3$

text: The question: > Find a 99% confidence interval for the difference in mean oxygen consumption > of an algae that is in 100% concentrated seawater with algae that is not in > 100% seawater (i.e., 50% or 75%). Is this the correct approach? > If we let $X_1$ be the random variable corresponding to the mean oxygen > consumption for algae in 100% seawater, $X_2$ for 50% and $X_3$ for 75% then > I want to find the standard error of $X_1+\frac{X_2+X_3}{2}$. If so, I know how to proceed from there. Otherwise, I don't know how to start this.


title: What if the MVUE depends on the parameter?

text: The minimum variance, unbiased estimator $\hat \theta$ of $\theta$ is defined by $$\hat \theta = \text{argmin}_{\hat \theta} \; \mathbb{E} \left( (\hat \theta - \theta)^2 \, | \, \theta\right), \quad \text{ subject to } \quad \mathbb{E} \left( \hat \theta \, | \, \theta\right) - \theta = 0 $$ where the expectation, $\mathbb{E}$, is taken over the data, $Y$ (given the parameter $\theta$). This yields an expression for $\hat \theta$ in terms of $y$. I.e. $\hat \theta = \hat \theta (y)$. But if instead we obtain $\hat \theta = \hat \theta (y; \theta)$, what do we do then? As far as I can tell, there's nothing in the definition of MVUEs above that guarantees independence of $\hat \theta$ from $\theta$. One could achieve such independence by letting the expectation be over $\theta$ as well as the data $Y$. But this assumes that $\theta$ is a random variable, and then we're doing Bayesian inference. What if we have a strictly frequentist perspective? For example, the ordinary least squares estimator, $(X'X)^{-1} X' \boldsymbol{y}$, is the BLUE (so the MVUE among linear, unbiased estimators) for the model $ \boldsymbol{Y} = X \boldsymbol{\theta} + \boldsymbol{\eta}$, and it does not have explicit dependence on $\theta$. Jolly good. Unfortunately I can't think of a case where the expression for the MVUE does depend on the parameter, however, again, there is nothing apparent in the definition of MVUEs above that can guarantee this. ## Edit If the requirement "for all $\theta$" is added the the definition of MVUEs (also called _uniform_ MVUEs (UMVUEs), that implies that there is not dependence on the parameter itself in the expression of the estimator. However, UMVUEs may not exist. As I mention in the comments, it seems weird that a frequentist would completely have to abandon the MVUE criterion in these "degenerate" cases, while a Bayesian could easily make sense of them too by taking expectations over the parameter as well, instead of requiring uniform optimality.


title: Is the product of two independent random variables independent of a random variable which is not independent of one of the factors?

text: Let $A$, $B$, $C$ be random variables. $A$ and $B$ are not independent $A$ and $C$ are independent. $B$ and $C$ are independent. Are $A$ and the product $BC$ independent? This is, is $E[ABC]=E[A]E[BC]=E[A]E[B]E[C]$??


title: Relation between observed power and p-value?

text: I am trying to understand the relation between observed power and p-value in Stephane's reply, which I think is based on J. M. Hoenig and D. M. Heisey (2001) "The Abuse of Power: The Abuse of Power Calculations for Data Analysis," The American Statistician 55(1), 19-24, ![enter image description here](http://i.stack.imgur.com/ByDeq.png) I don't quite understand the quote. I hope that someone could state it in a more clear way. Following are what I don't understand.   1. In the Z test, I was wondering why the parameter value $\delta$ is $\delta = \sqrt{n} \mu/\sigma$, but I think a parameter of a distribution shouldn't depend on the sample size $n$?   2. The $100(1-p)$th percentile $Z_p$ of the standard normal distribution doesn't depend on the sample $X$, does it? Why is it the observed statistic? Does the observed statistic mean the test statistic of the Z-test?   3. Why is the cdf $G_\delta$ of the p-value with the parameter value $\delta$ given by $$G_{\delta}(p) = 1 - \Phi(Z_p-\delta)?$$   4. The observed power depends on an observation $X$ of the sample, doesn't it? It is defined as the value of the power function at some estimate of the parameter $\delta$ from the sample $X$? Why is the observed power given by $$G_{Z_p}(\alpha) = 1- \Phi(Z_{\alpha} - Z_p)?$$   5. Does the relation between the observed power and the p-value for Z-test still apply for other tests? Thanks!


title: Are there methods to learn a projection method into euclidean space, given a set of pairwise distances?

text: Basically, right now I am trying to do some nearest-neighbour searching on an approach where I don't have points in Euclidean space. To do a nearest- neighbour search currently, I take the query and compare it (with a special similarity function) to each model in the database to find the most similar one. With this similarity function, I can compute pairwise distances between all points in my database, and I'm interested in whether it's possible to learn a projection into Euclidean space such that any new query can just be projected and nearest neighbour searches can then be performed in that Euclidean space. Does such a method exist, or am I hoping for a holy grail that hasn't been solved?


title: Run a Monte Carlo simulation of a startup while conceiving it? Has it been done?

text: As an entrepreneur you try to find a model of your business. If you simulate this via Monte Carlo you could see sensitive areas in e.g. the marketing strategy and how it affects e.g. sales. I like to think of it as a prediction device. My intuition is that a simulation allows for "what if" questions in a way shorter time. A strategy change that materializes to consequences in months can be examined in days or hours. Also in real life there is no "undo" button. My second intuition is that a startup is small enough to allow for gradually building the simulation in parallel to the real business. The validity of the simulation can be easly examined with the real world data. Literally simulate your startup. Has it been done before?


title: Accurate comparison over time period

text: I have 52 samples daily, weekly and biweekly data for 3 subjects over the course of a 10 month period. At two points during the 10 month period the resolution increases to daily samples as the people begin taking a supplement. I have been considering a percentage of supplement that is (total supplication/total all samples) but I see that that was flawed because the time of the total for all samples is so grossly greater than the time during the supplement. Now I am considering (total supplement/total before) and (average supplement/average before). I think that the average of the supplements and before eliminates the time factor problem. Here is a link to the sample dates: https://docs.google.com/spreadsheet/ccc?key=0Agl7ZLsLCsDQdFdzRGdhdU9Tb0hGNWtBOGkyMURnZlE Would someone help me make an accurate comparison?


title: What regression analysis should I perform on my data and why?

text: I am a law student researching which factors influence the CSR (corporate social responsibility, `GSE_RAW`) behavior of companies. As my studies didn't offer any statistics courses, I'm having trouble to understand what type of statistical analysis I should perform on my data. After describing the data, I hope some of you can tell me more about this. Two groups of possible factors / variables influencing CSR have been identified: company-specific and country-specific. First, company-specific variables are   * `MKT_AVG_LN`: the marketvalue of the company   * `SIGN`: the number of CSR treaties the company has signed   * `INCID`: the number of reported CSR incidents the company has been involved in Second, each of the 4,000 companies in the dataset is headquartered in one of 35 countries. For each country, I have gathered some country-specific data, among others:   * `LAW_FAM`: the legal family the countries' legal system stems from (either French, English, Scandinavian, or German)   * `LAW_SR`: relative protection the countries' company law gives to shareholders (for instance, in case of company default)   * `LAW_LE`: the relative effectiveness of the countries' legal system (higher value means more effective, thus for instance less corrupted)   * `COM_CLA`: a measurement for the intensity of internal market competition   * `GCI_505`: mesurement for the quality of primary education   * `GCI_701`: measurement for the quality of secondary education   * `HOF_PDI`: power distance (higher value means more hierarchical society)   * `HOF_LTO`: country time orientation (higher means more long-term orientation)   * `DEP_AVG`: the countries' GDP per capita   * `CON_AVG`: the countries' average inflation over the 2008-2010 timeframe In order to make an analysis on this data, I "raised" the country-level data to the company-level. For instance, if Belgium has a `COM_CLA` value of 23, then all Belgian companies in the dataset have their `COM_CLA` value set to 23. The variable `LAW_FAM` is split up into 4 dummy variables (`LAW_FRA`, `LAW_SCA`, `LAW_ENG`, `LAW_GER`), giving each company a 1 for one of these dummies. This all results in a dataset like this:               COMPANY MKT_AVG_LN ... INCID ... LAW_FRA LAW_SCA ... LAW_SR LAW_LE COM_CLA ... etc     ----------------------------------------------------------------------------------        1      1.54          55          0       1          34     65     53        2      1.44          16          0       1          34     65     53        3      0.11           2          0       1          34     65     53        4      0.38          12          1       0          18     40     27        5      1.98         114          1       0          18     40     27        .       .             .          .       .           .      .      .        .       .             .          .       .           .      .      .      4,000    0.87           9          0       1           5     14     18      Here, companies 1 to 3 are from the same country A, and 4 and 5 from country B. My DV, `GSE_RAW` is a numerical value for each companies' CSR behavior given by a rating agency.   * I believe the country-level variables are also called "categorical" variables, as many companies share the same value for these variables (in the example above, companies 1 to 3 all share the same values for `LAW_FRA` to `COM_CLA`). I believe to have found out that "categorical" variables are also known as fixed factors. Is all this true?   * I believe an OLS regression analysis is not the proper model here because of the categorical (country-level) variables. It has been proposed to use "Generalized Linear Models" (GLS), using the country-level variables as (fixed?) "factors" and the company-level variables as "covariates". Is this correct? And as a subquestion: why exactly is OLS not appropriate because of the country-level variables? What is it what they do in the OLS calculations that makes them set off the regression? `[edit 1]` I am using SPSS for statistical analysis `[edit 2]` Here my attempt to create a GLM using this data. However, I am unable to not get the "you haven't specified a custom model" Do I have to select all 4 variables here (becaus I want a beta and significance level for all 4 of them to construct a regression model)? And if so, why do I have to do this twice? I already said in a previous dialogue box that `DEP_AVG` and `CON_AVG` are fixed factors and that `SIGN` and `INCID` are covariates. Why would I, for instance, insert `INCID` here as a covariate, but not include it in the model building dialogue? Also, I really don't understand the output I'm getting, since it is very different from ordinary OLS output (the only output I'm slightly comfortable with).   * Am I now doing the right analysis?   * How can I get a regression model from this? ![enter image description here](http://i.stack.imgur.com/ewFeo.png)