title: What does this residuals versus fitted plot mean about my model?

text: I have a model that attempts to predict a nation's quality of life index by it's moral indifference to contraception and moral rejection of gambling. Initially the model contained several predictors, but I eliminated most using backwards elimination via AIC. Here is a summary of the model (generated using R):               > summary(fit1)          Call:     lm(formula = Quality.of.life.index ~ Morally.unacceptable.ga +          Not.a.moral.issue.co, data = qli_and_moral_ind)          Residuals:         Min      1Q  Median      3Q     Max      -89.670 -25.443  -4.732  36.129  64.441           Coefficients:                             Estimate Std. Error t value Pr(>|t|)         (Intercept)             143.1410    32.7499   4.371  0.00019 ***     Morally.unacceptable.ga  -1.7690     0.3603  -4.910 4.71e-05 ***     Not.a.moral.issue.co      1.4471     0.7925   1.826  0.07981 .       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 40.39 on 25 degrees of freedom     Multiple R-squared:  0.6079,    Adjusted R-squared:  0.5765      F-statistic: 19.38 on 2 and 25 DF,  p-value: 8.266e-06      There are two plots of the model that I can't interpret: ![Norm QQ Plot](http://i.stack.imgur.com/MBXmd.png)![residuals versus fitted plot](http://i.stack.imgur.com/NzRor.png) According to the web, the residuals plot above may indicate predictable error, ie that I'm missing some variable in my model. Is that assessment correct? If so, what should I consider adding to the model? It kind of looks like the $y = x^3 - x$ graph - maybe add a cubed term?


title: What does having "constant variance" in a linear regression model mean?

text: What does having "constant variance" in the error term means. As I see it , we have a data with variable and 1 independent variable. This is one assumption of linear regression. I am wondering what this homoscedasticity means. Since if I am having 500 rows I would have single variance value which is obviously constant. With what variable should I compare the variance?


title: Time series: Testing structural change of beta parameters of an OLS over time

text: I have a linear AR(2) regression model of an time series that looks like this: $y_{t} = \beta_0 + \beta_1 y_{t-1} + \beta_2 y_{t-2} + \beta_3 x_t + \epsilon_t$ I am modelling the contribution of the time series $\\{ y_{t-1}, y_{t-2}, x_t\\}$ on tomorrow's $y$, $y_{t}$. I want to test whether these contributions change over time. I do a rolling regression over a period of ~10 years giving me in total a set of about 2000 beta parameters $ \\{\hat{\beta_1}, \hat{\beta_2}, \hat{\beta_3} \\}$ over time. Indeed inspecting the change qualitatively reveals that there is a trend in the value the parameters over time. My question: How can I test formally that there is indeed a structural change in the parameters over time?


title: Vapniks proof of the basic lemma

text: In his book _Statistical Learning Theory_ (1998), Vladimir Vapnik proves an inequality needed to prove a bound on the risk for indicator loss functions. Theorem 4.1 on page 133 he derives the following inequalities $$I \ge \int_{Z_2} [\nu(\alpha^*,Z_2)-P(\alpha^*)-1/l]dF(Z_2)=P \left[\nu(\alpha^*,Z_2)\ge P(\alpha^*)^{-1/2} \right]$$ $$=\sum C_l^kP(\alpha^*)^k[1-P(\alpha^*)]^{l-k}$$ I absolutely don't get these two equalities. Furthermore, I think in the first integral there should be a $+ 1 / l$ Is anybody able to give me some idea? ++++++++++++++++++++++++++++++++++++++ EDIT: I NOW HAVE A REFORMULATION OF THIS PROOF MADE BY VAPNIK IN MY OWN WORDS AND NOTATION. THE LAST INEQUALITIES ARE MY PROBLEM Let $Z(2l)$ be the spaces of random selections of magnitude 2l. A element of this space we call $Z^{2l}=\left(z_1,\ldots,z_l,z_{l+1},\ldots,z_{2l}\right)$.\ For a loss function $\\{ L(z,\alpha) \\} _{\alpha \in \Lambda}$ we calculate the empirical risk $R_{emp}^{(1)} (\alpha) = \frac{1}{l} \sum_{i=1}^{l}L(z_i,\alpha)$ on $Z_1 = \left(z_1,\ldots,z_l \right)$ and the empirical risk $R_{emp}^{(2)} (\alpha) = \frac{1}{l} \sum_{i=l+1}^{2l}L(z_i,\alpha)$ on $Z_2 = \left(z_{l+1},\ldots,z_{2l} \right)$ von $Z^{2l}$. We define the random variables: $$\rho^\Lambda(\alpha, Z^{2l}) := | R_{emp}^{(1)}(\alpha)-R_{emp}^{(2)}(\alpha) |$$ $$\rho^\Lambda(Z^{2l}):= \sup \limits_{\alpha \in \Lambda} \rho^\Lambda (\alpha,Z^{2l})$$ $$\pi^\Lambda(\alpha,Z_1):= |R(\alpha)-R_{emp}^{(1)}(\alpha)|$$ $$\pi^\Lambda (Z_1):= \sup\limits_{\alpha\in\Lambda} \pi^\Lambda (\alpha,Z_2) $$ and assume $\pi^\Lambda (Z_1)$ and $\rho^\Lambda (Z^{2l})$ being measurable with respect to the probability measure $P$ generating the data. We want to show: $$P \left \\{ \pi^\Lambda (Z_1) > \epsilon \right\\} < 2 P\left\\{ \rho^\Lambda (Z^{2l}) > \epsilon - \frac{1}{l} \right\\}$$ By definition there exists $\epsilon_*:=\epsilon-\frac{1}{l}$ such that $$P \left \\{ \rho^\Lambda(Z^{2l}) > \epsilon_* \right\\}=\int_{Z(2l)} 1\left[ \rho^\Lambda (Z^{2l}) >\epsilon_* \right] dF(Z^{2l})\\\=\int_{Z_1} dF(Z_1) \int_{Z_2} 1\left[\rho^\Lambda (Z^{2l}) > \epsilon_*\right]dF(Z_2)$$ Let $Q$ be the event: $$ Q:= \left\\{ Z_1 : \pi^\Lambda (Z_1) > \epsilon \right\\}$$ We get $$ P \left \\{ \rho^\Lambda(Z^{2l}) > \epsilon_* \right\\} \geq \int_Q dF(Z_1)\underbrace{ \int_{Z_2(l)} 1 \left[ \rho^\Lambda (Z^{2l}) > \epsilon* \right] dF(Z_2) } _{Int}$$ We now examine the integral $Int$. Here $Z_1$ is fixed with $\pi^\Lambda (Z_1) > \epsilon$. That is why there exists $\alpha^* \in \Lambda$ such that $$ |P(\alpha^*) - R_{emp}^{(1)}(\alpha^*)|>\epsilon $$ for $$P(\alpha^*)=\int L(z,\alpha^*) dF(z)$$ with $$L(z,\alpha)=\frac{1}{2}|y-f(x,\alpha)|$$ then is $$Int = \int_{Z_2} 1\left[ \sup\limits_{\alpha\in\Lambda} \rho(\alpha, Z^{2l}) > \epsilon_* \right] dF(Z_2)\\\ \geq \int_{Z_2} 1\left[ \rho(\alpha^*, Z^{2l}) > \epsilon_* \right] dF(Z_2) $$ Without loss of generality we assume now $R_{emp}^{(1)}(\alpha^*)<P(\alpha^*)-\epsilon$ such that $$R_{emp}^{(2)}(\alpha^*) \geq P(\alpha^*)-\frac{1}{l}$$ is sufficient for $$\rho(\alpha^*, Z^{2l})=| R_{emp}^{(1)}(\alpha^*)-R_{emp}^{(2)}(\alpha^*)| > \epsilon_* $$ and we get $$ Int \geq \int_{Z_2} 1 \left[ \ R_{emp}^{(2)}(\alpha^*) \geq P(\alpha^*)-\frac{1}{l} \right]dF(Z_2)\\\ = P \left\\{ R_{emp}^{(2)}(\alpha^*) \geq P^{-1/2}(\alpha^*) \right\\}\\\ = \sum_{\frac{k}{l} \geq P(\alpha^*) - \frac{1}{l}} \binom{l}{k} \left[ P(\alpha^*) \right]^k \left[ 1-P(\alpha^*)\right]^{l-k}\\\ > \frac{1}{2}$$ resubstituting yields \begin{align*} P \left \\{ \rho^\Lambda(Z^{2l}) > \epsilon_* \right\\} > \frac{1}{2} \int_Q dF(Z_1)= \frac{1}{2} P\left\\{ \pi^\Lambda (Z_1) > \epsilon \right\\} \end{align*} qed


title: What are the 'big problems' in statistics?

text: Mathematics has its famous Millennium Problems (and, historically, Hilbert's 23), questions that helped to shape the direction of the field. I have little idea, though, what the Riemann Hypotheses and P vs. NP's of statistics would be. So, what are the overarching open questions in statistics? **Edited to add:** As an example of the general spirit (if not quite specificity) of answer I'm looking for, I found a "Hilbert's 23"-inspired lecture by David Donoho at a "Math Challenges of the 21st Century" conference: High-Dimensional Data Analysis: The Curses and Blessings of Dimensionality So a potential answer could talk about big data and why it's important, the types of statistical challenges high-dimensional data poses, and methods that need to be developed or questions that need to be answered in order to help solve the problem.


title: Variable selection for regression - the subselect package

text: No regular here will be unaware of the perils of using stepwise and similar automatic methods for variable selection in regression analysis. But preferred alternatives, such as the lasso or elasticnet, have there own difficulties. I can't find anywhere in the archive here a discussion of the methods provided by the **subselect** package in R, which I've just come across - the package has, so far as I can see, existed for a decade and more, and presumably has proved useful. In addition to a variation on the leaps procedure, **subselect** offers three algorithms (which it calls anneal, genetic and improve) for variable selection for different kinds of analysis. Have these procedures (or any of them) proved to be of value in variable selection?


title: Factor analysis of questionnaires composed of Likert items

text: I used to analyse items from a psychometric point of view. But now I am trying to analyse other types of questions on motivation and other topics. These questions are all on Likert scales. My initial thought was to use factor analysis, because the questions are hypothesised to reflect some underlying dimensions.   * But is factor analysis appropriate?   * Is it necessary to validate each question regarding its dimensionality ?    * Is there a problem with performing factor analysis on likert items?   * Are there any good papers and methods on how to conduct factor analysis on Likert and other categorical items?


title: How to get Pareto IV parameter estimates

text: I have a serie with 850 observations, and I need to fit the Pareto IV distribution. How could I do this in R? I read the guide VGAM, however, I'm not able to run it. If anyone knows, please provide me with a detailed (step by step) answer so that I can understand it...


title: A transform to change skew without affecting kurtosis?

text: I am curious if there is a transform which alters the skew of a random variable without affecting the kurtosis. This would be analogous to how an affine transform of a RV affects the mean and variance, but not the skew and kurtosis (partly because the skew and kurtosis are defined to be invariant to changes in scale). Is this a known problem?


title: Is there a biglm equivalent for coxph?

text: I am trying fit a Cox Proportional Hazard model with coxph with a cmoputer of 12GB of RAM. It keeps running out of memory. Is there a biglm equivalent of coxph?


title: Cramér's $V$ on Rao-Scott adjusted Pearson $\chi^2$

text: I have survey data with design weights for stratified sampling. My ultimate goal is to estimate Cramér's $V$ for contingency tables, a Pearson $\chi^2$ based measure. To account for the weights, I am thinking about using Rao-Scott adjustment to estimate $\chi^2$ statistics. However, then it is not clear how to get to Cramér's $V$. Is it accurate to estimate $V$ from a first-order Rao-Scott adjusted $\chi^2$, just like I would estimate it from an unadjusted $\chi^2$? Or should $V$ be adjusted as well? Thank you! Some notes: http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.aos/1176346391 (Original Article) http://www.amstat.org/sections/srms/proceedings/y2007/Files/JSM2007-000874.pdf (Some summary)


title: How to link the SE(p) of the binomial distribution $B(n,p)$ to the SE(coefficient) in GLM?

text: This is a general question regarding GLM with binomial distribution. I use the following data (with $N=400$ observations) as an example (predict toxicity using treatment):               mod1 <- glm(toxicity~treatment, family=binomial, data=Dat)      The description of the model can be written as: $Y_{i}=Binomial(1,p_{i})$ $logit(p_{i})=\alpha+\beta \times treatment_{i}$ where $y_{i}$ is 1 if the $i$th observation gives toxicity. The model output shows:                                 Estimate    Std. Error z value   Pr(>|z|)         (Intercept)        -0.140784   0.259487  -0.543    0.58744         treatment_A        -1.088235   0.254966  -4.268    1.97e-05 ***      I know that the SE of the coefficients are defined as the diagonal of the covariance matrix of the coefficients, while the covariance matrix is the inverse of the negative of the second partial derivatives of the Log- likelihood with respect to the coefficients, evaluated at the values of the coefficients that maximize the likelihood. On the other hand, from the binomial distribution, the SE of $p_{i}$ is $sqrt(p_{i}\times(1-p_{i})/n)$ ($n$ refers to the number of observations at $treatment_{i}$). According to the model, the sampling error of the coefficients (SE of coefficient) should come from the sampling error of the binomial distribution. My question is: 1) Is there any relations, and if so, how is the relation between SE of the coefficients and the SE of $p$ from the binomial distribution? 2) This actually is the motivation of my question. Looking at the binomial distribution only, and if all the observations are i.i.d., increasing the number of observations will improve the precision of {p}. And I can imagine the corresponding SE(coeff) will be decreased, thus increasing the power of the test. And that's why I would like to know how these are related.


title: Calculating statistical power

text: As I understand it, I need to know at least three aspects (out of four) of my proposed study in order to conduct power analysis, namely:   * type of test - I intend to use Pearson's r and ANCOVA/Regression - GLM   * significance level (alpha) - I intend to use 0.05   * expected effect size - I intend to use a medium effect size (0.5)   * sample size Could anyone recommend a good online power calculator that I can use to do _a priori_ power calculation. (Can SPSS do _a priori_ power calculation?) I have come across GPower but I am looking for a simplier tool!


title: Do you think smoothing constant value, alpha, in SES method is a control parameter or process parameter?

text: There is a debate in selecting the smoothign constant in Single Exponentioan Smoothing method by practitioner or considering it as a process parameter? Could you please provide your opinion regarding this issue?


title: How to use PLINK for repeated measures?

text: I am trying to use PLINK to do a GWAS for a outcome with repeated measures, namely each subject has time1 and time2 observations. Does anybody have some experiences or suggestions for this. I am thinking use a gxe model such as:               plink --bfile mydata --pheno outcomes.txt --gxe --covar time.txt --out results      If so, each subject in the outcomes.txt will have two measures, such that:               FID IID Pheno     1001 11  0.2     1001 11  0.4     1002 12 0.6     1002 12 0.8     ...      and use a time.txt like:               FID IID time     1001 11  1     1001 11  2     1002 12 1     1002 12 2     ...      I am wondering whether it is okay to do this way.


title: Dummy interaction variables are always non-stationary?

text: I want to know why we can include dummy interaction terms into time series models if they're always non-stationary? For example let $X_t$ be $I(0)$, $X_t \sim N(\mu,\sigma^2)$ and $D_t \in \\{0,1\\}$. Suppose the estimation window is $200$ and $D_t =0$ for $t=1,...,100$ and $D_t =1$ for $t=101,...,200$. The interaction is $D_tX_t$. We know that   * $E[D_tX_t|t\in\\{1,...,100\\}]=0$   * $E[D_tX_t|t\in\\{101,...,200\\}]=\mu$ Therefore it's non-stationary? The variance also changes for this variable at $t=101$. This of course holds for intercept dummy variables ($D_t$ by themselves); however in this case the variance is stationary (but the mean changes).


title: Replicating integration results from a paper

text: I’m reading Wacek’s paper _Parameter Uncertainty in Loss Ratio Distributions and its Implications_ and trying to figure out how to replicate some of the results. Table 6, on page 190 of the paper, contains some results. Under Model B1, the density, $f_x(x|\theta)$, is assumed to be Normal with known $\mu = 0.7067$ and $\sigma=0.0745$. The formula of density is given by Eqn (2.3). The pure premium is given by Eqn (5.1): $$\int_R^{L+R} \\! (x-R)f_x(x) \, \mathrm{d} x + L\int_{L+R}^\infty \\! f_x(x) \, \mathrm{d} x$$ where L (Limit) = 0.5, and R (Retention) = 0.7. According to the results provided in Table 6, after integrating, I'm supposed to get 2.02%, but what I'm getting is different. Here is the R code I'm using:               integrand1 <- function(x)       {(x-0.7)*(1/(0.0745*sqrt(2*pi)))*exp(-0.5*((x-0.7067)/0.0745)^2)}     l1 <- integrate(integrand1, lower = 0.7, upper = 0.75)          integrand2 <- function(x)       {0.05*(1/(0.0745*sqrt(2*pi)))*exp(-0.5*((x-0.7067)/0.0745)^2)}     l2 <- integrate(integrand2, lower = 0.75, upper = Inf)      Any thoughts?


title: Error amplification in linear regression with uncertainties in design matrix

text: I am running a least squares regression, and I am trying to estimate the uncertainty on the fit solution. The problem is, I have an uncertainty in my design matrix coefficients as well as in my predictors. So given $||Ax - b||_2$ being minimized, and uncertainties $\sigma_b$ and $\sigma_A$ on the coefficients of $b$ and $A$, what is the uncertainty $\sigma_x$ on the resulting vector $x$?


title: Analysis of binary variables

text: I have a data set consisting of about a quarter-million objects, each of which may have any of 30 particular features. So I might have Object 1: feature 3, feature 7   Object 2: feature 3, feature 29, feature 30   Object 3: feature 3, feature 7   Object 4: feature 1, feature 18, feature 20, feature 28   ... At this point I'm just doing exploratory analysis. I'm interested in seeing how the different features relate: does 29 always appear with 30? Does feature 7 often occur with feature 3? Can feature 10 be predicted from features 1, 2, and 3? Etc. What types of analysis are appropriate here? The problem is different from others I have worked on in the past because the state space is both small (each object 'holds' merely 30 bits of easily-compressible information) but paradoxically also large ($30 \choose 2$ is already too large to readily present, so even pairwise interactions are hard to study, and $30 \choose 3$ is much too large). My immediate thought is something that would identify which pairs, triples, quadruples, etc. have 'interesting' interactions, even though examining them individually is infeasible. But perhaps there are more interesting things to do? Basic references may be appropriate here (and would be appreciated if relevant).


title: Naive Bayesian Algorithm in R/SAS for categorical input variables?

text: Could anyone please let me know how to implement Naive Bayesian Algorithm in R or SAS?I have got a training dataset with all the categorical predictors and target variable(3 levels).I got to build a model and apply it on a different test dataset along with the probable target and its predicted probability. To be more clear,my first dataset 'A' contact 4 categorical input variables a,b,c,d and a target class 'T' of 3 levels.I need to train the model for this dataset initially.Then,I have got one more dataset 'B' with input categorical variables w,x,y,z and I need to predict the probable target class 'S' along with its probability here based on my previous built model.I want the entire thing to be done in R or SAS but couldn't find much resources.Sorry,if the question has been repeated.


title: Compound poisson process: Average size of claim will exceed £110

text: "An insurance company receives claims at a rate of two per week, the size of a claim in pounds having mean 100 and standard deviation 50. Assuming the compound Poisson process as a model, and using normal approximations, find the probabilities that, in a period of 50 weeks the average size of claim will exceed £110". As it's a PP, I have $$Q(t) = \sum_i^{N(t)} Y_i$$ I also have $\mu = 100$ and $\sigma = 50$. I want to work out $P\left(\frac{Q(t)}{N(t)} > 110 \right) \implies P\\{Q(t) - 110N(t) > 0\\}$ This tells me I have a new compound PP which can be written as $$ \tilde{Q(t)} = \sum_i^{N(t)} (Y_i - 110) $$ Now we want $E[\tilde{Q(t)}]$ and $Var[\tilde{Q(t)}]$ $$ E[\tilde{Q(t)}] = E[Q(t) - 110N(t)] = E[Q(t)] - 110E[N(t)] = \lambda t \mu - 110 \mu = -1000$$ $$ Var[\tilde{Q(t)}] = Var[Q(t)] - 110Var[N(t)] = \lambda t (\sigma^2 + \mu^2) - 110\sigma^2$$ However, this variance is wrong. What is the correct way to work out the variance?


title: Imputation with R and MICE

text: A short `mice`-related question as follows: when running simple mice imputation the function goes through well. But `complete()` function does not substitute my `NA`'s at all. I guess it must be something wrong with my dataset but I can't see what. It is just a matrix as any other although with roughly equal number of rows and columns. Maybe this is the source of the problem. What should I look into?


title: Parametric Surface Reconstruction from Contours with Quick Rescaling

text: I'm looking to construct a 3-D surface of a part of the brain based on 2-D contours from cross-sectional slices from multiple angles. Once I get this shape, I want to "fit" it to another set of contours via rescaling. I'm aspiring to do this in the context of an MCMC analysis (So as to be able to make inferences, so it would be very nice if I could easily compute the volume of the rescaled surface, and the minimum distance between a given point of the distance. (Accurate approximations are fine). What would be a good image reconstruction algorithm that allows for volume and distance to be quickly calculated?


title: Is my weatherman accurate?

text: A question which bothered me for some time, which I don't know how to address: Every day, my weatherman gives a percentage chance of rain (let's assume its calculated to 9000 digits and he has never repeated a number). Every subsequent day, it either rains or does not rain. I have years of data - pct chance vs rain or not. _Given this weatherman's history_ , if he says tonight that tomorrow's chance of rain is X, then what's my best guess as to what the chance of rain really is?


title: Based on my data, is Jack likely to be clumsy?

text: I posted this on mathoverflow, but they sent me here. This question relates to a problem I had at work a while ago, doing a little data mining at a car rental company. Names changed, of course. I'm using Oracle DBMS if it matters. There was a flight of steps out the front of our building. It had a dodgy step on it, on which people often stub their toes. I had records for everyone who works in the building, detailing how many times they climbed these steps and how many of these times they stubbed their toes on the dodgy step. There's a total of 3000 stair-climbing incidents and 1000 toe-stubbing incidents. Jack climbed the steps 15 times and stubbed his toes 7 times, which is 2 more than you'd expect. What's the probability that this is just random, vs the probability that Joe is actually clumsy? I'm pretty sure from half-remembered statistics 1 that its something to do with chi-squared, but beats me where to go from there. ... Of course, we actually had several flights of steps, each with different rates of toe stubbing and instep bashing. How would I combine the stats from those to get a more accurate better likelihood of Joe being clumsy? We can assume that there's no systematic bias in respect of more clumsy people being inclined to use certain flights of steps.


title: How to learn how to use a new statistical GUI?

text: The question is in the header, but I would extend the context a bit. Next semester I am due to be a TA in a course in statistics, where I would need to help sociology students learn to use SPSS. I don't know SPSS, yet, and would like to learn how to use it. I was thinking of taking a simple dataset, and start reviewing it with methods I know, thus starting to map out where are methods I know of. And once finished, to try and explore more options. Can someone propose other/better strategies to master a new statistical GUI? (in my case SPSS, but it could apply to many other GUI's)


title: Choosing the appropriate method to determine risk factor (logistic regression)

text: My question is about logistic regression, and I want you to advise me to use the appropriate method for my problem. Here is the description: My goal is to determine the risk factors for a disease (a virus) which has many types. One individual can be infected with many types, and every type has a "rated risk" (high risk: `HR`, intermediate risk: `IR`, low risk: `LR`). So I have a database with 325 individuals, and for every individual I have one column that indicates his state (`0`: not infected and `1`: infected), 3 other columns `HR`, `IR` and `LR` (`0`: false, `1`: true, for every column). I have also many other columns indicating the exact infections: `Inf1`, `Inf2`, `Inf3`,..., `Infk` (`0`: not infected, `1`: infected, for every column). As risk factors, I have many other social variables (columns), which can be quantitative or qualitative... One of them is the Sex activities in the last 12 month (0:No ; 1:Yes). If the response is Yes, then I have the exact number of partners. How can I enter these two variables in the model, as they are related to each other. My sample contains only 13% (43/325) infected individuals.   Is it considered a rare event? Will my estimate be biased?   My question: what is the best method to use to determine the risk factors? * * * I have done a binary logistic regression (using the state as dependent variable) and the model was not able to predict the positive individuals.


title: Calculating scale mean from item means

text: I'm working on a review paper and need to collect the means and standard deviations of a given measure (such as a measure of depression) from papers of interest. However, some authors report means and standard deviations for each item on the measure, but do not calculate the overall mean and standard deviation. For example, for a measure with a 1-5 Likert-type response scale, rather than reporting a mean depression score of 2.5, sd = .34, for the whole scale, they'll report: item 1 mean 2.4, sd=.41, item 2 mean 2.5, sd=.38, etc. 1) Can I calculate the scale mean by simply taking the average of the item means? 2) Is it possible to calculate the accompanying scale standard deviation? And if so, how? It's probably not as simple as taking the average of the item standard deviations.... If it isn't possible, I can just contact the authors, but I hate to bother them if the needed stats can be calculated from the provided information. Any advice would be greatly appreciated! Thank you!


title: Classification measures for linear classifier

text: Let $\mathcal{H}\colon\mathbf{w}\cdot\mathbf{x}+b=0$ be a separating hyperplane, which some binary linear classifier results in. Let $\mathbf{x}_t$ be an unseen, new sample that appears and needs to be classified. We can predict the truth label of $\mathbf{x}_t$ simple by computing the sign of the distance between $\mathbf{x}_t$ and $\mathcal{H}$, i.e., $$ y_t=\operatorname{sgn}(d_t), $$ where $$ d_t=\frac{\mathbf{w}\cdot\mathbf{x}_t+b}{\lVert\mathbf{w}\rVert}. $$ Then, a degree of confidence (probabilistic estimation), $s_t$, about the predicted label may be computed using the sigmoid function $S\colon\mathbb{R}\to(0,1)$, given by $$ S(t)=\frac{1}{1+e^{-t}}. $$ That is, $s_t=S(d_t)$. May anyone give some alternative measure of $s_t$? I need to use $s_t$ as a probabilistic estimation for classification. Thanks in advance for any helpful comment.


title: Question about DIANA algorithm and plotting the result

text: I am using the diana() function in R from the "cluster" package to cluster some data. To plot the result I use pltree(): pltree(diana(data.dist)) So the problem is that I have 3 brances in one node (4th split). It seems confusing to me, because DIANA should divide cluster into two at each step. I tried to convert the result of using diana() into dendrogram object and to plot it with plot(): plot(as.dendrogram(daisy(data.dist))) Is it modification of DIANA algorithm in R or just an error of plotting? ![enter image description here](http://i.stack.imgur.com/QaAZX.png) Thanks!


title: P value and confidence interval for two sample test of proportions disagree

text: I'm using R to calculate the two-sample test for equality of proportions, where the two proportions are 350/400 and 25/25. So:               > prop.test(c(350,25),c(400,25))                                                                                2-sample test for equality of proportions with continuity correction          data:  c(350, 25) out of c(400, 25)      X-squared = 2.4399, df = 1, p-value = 0.1183     alternative hypothesis: two.sided      95 percent confidence interval:      -0.17865986 -0.07134014      sample estimates:     prop 1 prop 2       0.875  1.000           Warning message:     In prop.test(c(350, 25), c(400, 25), correct = FALSE) :       Chi-squared approximation may be incorrect      What I can't reconcile on my own is that the p-value is greater than 0.05, and yet the 95% confidence interval for the difference does not include 0. I thought there was an 'if and only if' relationship between the two (The p-value < alpha iff the (1-alpha) confidence interval of the difference does not include 0). What am I not seeing? My only guess is there's something fundamental that I'm misunderstanding, or that it has something to do with that warning message about chi-squared approximation.


title: Struggling with copula theory

text: I'm really struggling with bivariate copula's. Long story short, I can only use Gaussian copulas. I'm therefore interested in the joint PDF for which the Gaussian copula can be applied. So for example:   * The Gumbel copula is used for extreme distributions.   * The Gaussian copula is used for linear correlation.   * The Archimedean copula and the t-copula are used for dependence in the tail. So this clearly "limits" the joint PDFs for which the Gaussian Copula can be used to accurately model, and there are many cases that it simply won't apply. I.e., the joint PDF must have linear correlation, small or no tail dependence, and it is not ideal for the extreme distributions. Are there any more restrictions? Furthermore, does anyone know of a simple joint PDF that can be accurately used? I'd ideally like to find one which involves Uniform, Gamma and Beta distributions. Several sources I've looked at just plough through this and state that the Gaussian copula can be used all over the place–but I just can't see how?! Have these authors just gotten carried away and ignored the limitations?


title: use of weights in svyglm vs glm

text: I would like to know how the treatment of weights differs between `svyglm` and `glm` I am using the `twang` package in R to create propensity scores which are then used as weights, as follows (this code comes from the `twang` documentation):               library(twang)     library(survey)     set.seed(1)          data(lalonde)          ps.lalonde <- ps(treat ~ age + educ + black + hispan + nodegree + married + re74 + re75,      data = lalonde)          lalonde$w <- get.weights(ps.lalonde, stop.method="es.mean")     design.ps <- svydesign(ids=~1, weights=~w, data=lalonde)          glm1 <- svyglm(re78 ~ treat, design=design.ps)          summary(glm1)          ...     Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)   6685.2      374.4  17.853   <2e-16 ***     treat         -432.4      753.0  -0.574    0.566          Compare this to:               glm11 <- glm(re78 ~ treat, weights=w , data=lalonde)     summary(glm11)          Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)   6685.2      362.5  18.441   <2e-16 ***     treat         -432.4      586.1  -0.738    0.461        So the parameter estimates are the same but the standard errors for the treatment are quite different. How does the treatment of weights differ between `svyglm` and `glm` ?


title: Pattern recognition learning by examples

text: I am learning pattern recognition. But in the book I am reading there is almost exclusively a theory with few examples. Is there a special book / web site etc, with a lot of tasks covering all main topics of pattern recognition? I understand the theory, but I need practice so I'd like to find a book with list of problems like books on mathematical analysis and so on.


title: Computing the dyadic add portion of the Walsh-Fourier spectral density

text: I am attempting to build a Walsh-Fourier spectral density and it appears that it is first required to compute the logical covariance which in turn involves a dyadic add. I am not at all familiar with dyadics or their operations and introductory references are hard to come by. In Stoffer (1988) the logical covariance of a categorical series $X(0), X(1),\dots,X(N-1)$ is described as being: \begin{align} \tau(j)=N^{-1} \sum_{j=0}^{N-1} \gamma(j\oplus k-k) \end{align} where $j\oplus k$ is the dyadic addition. $\gamma$, is our usual autocovariance, $\gamma(h)=cov\\{X(n), X(n+h)\\}$. The Walsh-Fourier spectral density is then: \begin{align} f(\lambda)=\sum_{j=0}^{\infty}\tau(j)W(j, \lambda) \end{align} where $W(j, \lambda)$ is the $j$th sequency (zero-crossings) with $0\leq\lambda < 1$. I'm sure a HMM would be great for categorical times series but at the moment I am restricted to spectral analysis so I must continue with this approach. It's not exactly homework. It's a final project that has moved a little beyond the coursework. The professor is a little hard to get a hold of and thus the question is posed here. Is it just addition? I haven't even been able to confirm that.


title: Comparison between bootstrap percentiles confidence interval and BCa confidence interval?

text: I have generated a nonparametric percentiles bootstrap confidence interval ($32.27143, 51.08571$), and a BCA confidence interval ($33.26, 53.49$) with an initial sample size of $n=7$. Evidently, the BCA interval is larger than the percentiles interval. Shouldn't the BCA interval be better than the percentiles interval (and by "better" I mean "smaller")? I noticed also that when the sample size is larger, for example $n=40$, BCa confidence interval decreases.


title: Calculating level 2 variable based on mean of level 1 variable with multilevel model in R

text: I have 299 surveys collected from 299 individuals working at 26 different locations. I want to understand how the location specific features relate to the individual survey scores. The only inference I have as to location features is gathered from the individual survey scores. Is it a valid strategy to calculate means for each location based on the individual scores, and include this as a level 2 variable? Further, does it also make sense to include the same variable but as the level 1 variable, with slope varying freely between locations, if I want to compare the relative usefulness of the mean (best estimate of 'reality') to a persons individual score? (their perception of reality and response biases). I feel like I may have some circularity in the logic. My implementation in R for one of the variables of interest follows, any feedback is welcome!               lmer(X21~X25+meanX25+(X25|X1),data=datai)     Linear mixed model fit by REML      Formula: X21 ~ X25 + meanX25 + (X25 | X1)         Data: datai        AIC  BIC logLik deviance REMLdev      1079 1105 -532.7     1056    1065     Random effects:      Groups   Name        Variance Std.Dev. Corr         X1       (Intercept) 0.384983 0.62047                        X25         0.012382 0.11127  -1.000       Residual             1.936068 1.39143      Number of obs: 299, groups: X1, 26     Fixed effects:                 Estimate Std. Error t value     (Intercept)  1.13616    0.38013   2.989     X25          0.56683    0.05265  10.766     meanX25      0.33897    0.12213   2.775          Correlation of Fixed Effects:             (Intr) X25        X25     -0.119            meanX25 -0.838 -0.389


title: Regression tree model splitting too far - random data column?

text: I am attempting to make some regression trees with many potential independent variables which comprise both categorical and continuous data types of widely varying scales i have been using a few different methods of modelling with different validation methods ANNs, Decision trees, Boosted trees. Validation with hold back, 3 way hold back and k-fold i have decided to focus on decision trees due to the simple ability to convey results to others. I have also decided that due to holdback validation producing some quite different models dependent on which members are used for the validation that the dataset is too small for this What i am left with is k-fold decision tree which appears to split too far despite no k-fold validation drop off . I was shown a method (by a SAS rep) where a variable is inserted into the model making process that is essentially junk to help identify this using this method i have found that the 'junk' variable gets incorporated into the model sometimes quite early on and splits are performed beyond this. Is it valid to prune the tree back to remove these to be left with a more true model? incidentally, the hold back validation method produces similar models (in terms of number of splits) as this k-fold pruned model sorry for the long read, help would be appreciated! Edt: I would like to add that when using 'junk' data in holdback it never gets incorporated into the model due to the early stopping that occurs in splitting variables


title: Why am I not able to fit a zero inflated poisson distribution?

text: Following what is suggested here http://stackoverflow.com/questions/7157158/fitting-a-zero-inflated-poisson- distribution-in-r               > stat         x    N      1: 0  478      2: 1  901      3: 2 1101      4: 3  873      5: 4  583      6: 5  250      7: 6   97      8: 7   31      9: 8   10     10: 9    2          # vect <- rep(stat$x, stat$N)     count <- c(478, 901, 1101, 873, 583, 250, 97, 31, 10, 2)     vect <- rep(0:9, count)     library(fitdistrplus)     library(gamlss)          fit <- fitdist(vect, "ZIP", start=list(mu=2.4, sigma=0.1))     # mu = 2.64, sigma = -0.14, log = TRUE): sigma must be between 0 and 1       The plots are from regular poisson fit. As I see there are more zeroes, and gof is 0.00087, so I hope ZIP could help. ![enter image description here](http://i.stack.imgur.com/zDK5d.png) However, if I use `zeroinfl` from `pscl`               summary(zeroinfl(x ~ 1, dist="poisson", data=data.frame(x=vect))          Pearson residuals:         Min      1Q  Median      3Q     Max      -1.4945 -0.8607 -0.2269  0.4069  4.2096           Count model coefficients (poisson with log link):                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  0.88120    0.01134   77.73   <2e-16 ***          Zero-inflation model coefficients (binomial with logit link):                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  -3.7452     0.2597  -14.42   <2e-16 ***     ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1           Number of iterations in BFGS optimization: 10      Log-likelihood: -7853 on 2 Df          mu = exp(0.8812) = 2.41     zero = logit(-3.7452, inverse=T)=0.02308537


title: Bayesian justification for AIC/BIC

text: Can someone point me to a straightforward and comprehensible Bayesian discussion justifying AIC and/or BIC? Or even better, can someone give a self- contained such discussion in this forum?


title: Will two estimators converge to the same answer?

text: Say I have two estimators for the same quantity and using the same model, $E[f(X)]$. I also know that these two estimators are consistent, meaning, if we have a lot of data, they will be close to the true $E[f(X)]$ -- but all of this under the assumption that the model is _correct_! MLE is one kind of estimator like this. Will they converge to identical answers even if the model is not correct (i.e. the i.i.d. samples are not sampled from the parametric family of the model)?


title: Count data with one factor level containing only zeroes

text: I have a simple poisson glm with one predictor that has three levels. Unfortunately, for one level my response, the variable has only counts of zero. I expected very low counts (perhaps a one or a two at most). Because of one level having all zeros, neither my poisson glm, nor "zeroinfl" (from pscl), work. The glm throws out a massive standard error for the all zero factor level and the zeroinfl gives me an error message. Is there a statistical way to say that the other two groups are different from the one with all zeros?


title: Robust parameter estimation for shifted log normal distribution

text: I have a data set which fits a logNormal distribution quite well. (From a theoretical point of view, it is some hard-to-tackle quotient distribution). However, the data is quite dirty, so parameter estimation is far from trivial. Right now, my approach is the following:   1. Shift the distribution such that the minimum is almost 0.   2. Logspace the data   3. Use a robust Median and MAD parameter estimation (see Estimating parameters of a normal distribution: median instead of mean? for details) The result is significantly better than before (Maximum difference from empirical CDF 0.034 instead of 0.081 and 0.224 without using MAD). It's not perfect in particular on the long tail where I expect outliers. The additional location parameter helped a lot. However, using the minimum is a very crude heuristic. I obviously cannot expect to observe the true minimum, but depending on the sample size the observed minimum will always be some small x larger. Do you know any _robust_ parameter estimation method (+ a reference if possible) for the $e^{\mathcal{N}(\mu, \sigma)} + c$ distribution family? Note that e.g. scipy.stats.lognorm does also have such an additional third location parameter, just like the one I'm using, but I'm working in Java with my own code. **Update** : I've just come across a thesis on this topic:   * Estimating the Parameters of the Three-Parameter Lognormal Distribution   Rodrigo J. Aristizabal which includes pointers to some relevant literature, in particular to   * Estimating Parameters of Logarithmic-Normal Distributions by Maximum Likelihood   A. C. Cohen, Jr. but I find it hard to get a formula out of these publications that I could implement.


title: How do I conduct an analysis of an interrupted time series design with a comparison group?

text: I am not sure if this is the right place to ask this question (if not, please refer me to the right place), but it would be nice if someone could help me: **Description of the Data/the Scenario:** I have access to a huge amount of secondary data (>10000 cases) from a sourcing department of a company in the time-span from 2006-2012. The dataset consists of the ID of the part that has been sourced, the year of price-negotiation, the price resulting from the negotiation and the participating companies in the negotiation. In this data, I want to test, whether the participation of low-cost competitors (in price-negotiation) had an effect on the price-trend/level over the (following) years. Therefore, I plan to identify all parts that got negotiated (every year) from 2006-2012, then I want to split them into 2 groups, group 1 has low-cost country participation in negation after 2009 and group 2 has no low-cost country participation in negotiations at all (from 2006-2012). (Additional groups may be possible, but first I want to know how to do this with only two groups.) In essence, I have a longitudinal design with control group (respectively, I have a multiple group short interrupted time series quasi-experiment or Interrupted Time Series Design with Comparison Group) **Research questions:** Has the participation of low-cost countries (long-term) effects on the price levels? **The Design:**                Year:  06  07  08  09  10  10  11  12      Group1: O   O   O   O   X   O   O   O      Group2: O   O   O   O   -   O   O   O      **(Example) Dataset:**               ID   Time  Price  group      1    2006   1,5$  1     1    2007   1,5$  1     1    2008   1,5$  1     1    2009   1,3$  1     1    2010   1,2$  1     1    2011   1,2$  1     1    2012   1,1$  1     2    2006    10$  2     2    2007   9,9$  2     .    .      .     .      **Questions:** A) Which analyses-approaches do I need to perform (e.g.LGM, HLM, PS--> see Braver & Braham 2005)? B) Which programs (I have SPSS, otherwise preferable open source) do I need? C) How to code my data that they fit in the program/ fit the analysis-approach (long format or wide format, dummies)?


title: Covariate shift correction: standard implementation in R?

text: I have a dataset I'm working on that has some co-variate shift between the training set and the test set. I'm trying to build a predictive model to predict an outcome, using the training set. So far my best model is a random forest. How can I deal with the shifted distributions in the training vs. test set? I've come across 2 possible solutions that I've been able to implement myself:   1. Remove the shifted variables. This is sub-optimal, but helps prevent my model from over fitting the training set.   2. Use a logistic regression to predict whether a given observation is from the test set (after balancing the classes), predict "test set probabilities" for the training set, and then boostrap sample the training set, using the probabilities for sampling. Then fit the final model on the new training set. Both 1 and 2 are pretty easy to implement, but neither one satisfies me, as #1 omits variables that might be relevant, and #2 uses a logistic regression, when my final model is tree-based. Furthermore, #2 takes a few paragraphs of custom code, and I worry that my implementation may not be correct. What are the standard methods for dealing with covariate shift? Are there any packages in R (or another language) that implement these methods? /edit: It seems like "kernel mean matching" is another approach I could take. I've found lots of academic papers on the subject, but no one seems to have published any code. I'm going to try to implement this on my own, and will post the code as an answer to this question when I do.


title: Marginal posterior and prior are similar (and flat!)

text: I designed a Bayesian model and sampled the posterior using a MCMC algorithm. My problem is that the posterior marginal distribution of a given latent intermediate variable appears to be uniform just as the prior I assigned to it. In practice this variable is supposed to have a substential importance on the model. Moreover, the posteriors over the other variables sound in accordance with intuition and essentially monomodal. I am a bit confused with that situation. How to interpret this result? Have I to change my model? Considering that the problematic variable is an auxiliary variable that is not interpreted after inference, can I satisfied myself from these results, or can it be interpret as a failure of the modelisation ?


title: Normal distribution in R

text: I am trying to check what is the probability that a new observation is anomalous or not? Suppose I have the following set of observations:               x <- c(11,22,3,4,25,6,7,1,1,2,1)      where the sample mean of $x$ is 7.545455 and the sample standard deviation is 8.489566 Suppose I observe a new value of $x$, 111, which is anomalous because its value is considerably more than the sample mean of the normal distribution. How can I check in R that it is anomalous?


title: Sampling distribution of means

text: Im working through a textbook atm and Ive read that the mean of "the sampling distribution of means" is the same as the "populations" mean. To construct the sampling distribution of means, i need to take X * number of independent samples (where each sample is composed of >= ~30 observations) from the population, calculate the mean of each sample and plot the mean. Q1: Whats the minimum value that X (the number of means used to construct the sampling distribution) should be, or do you typically plot the distribution and increase X until it begins to look normal?


title: Visually plotting multi dimensional cluster data

text: I have a data set with 16 variables, and after clustering by kmeans, I wish to plot the two groups. What plots do you suggest to visually represent the two clusters?


title: How to test  a logistic regression model developed on a training sample on the data left out using R?

text: I have the summary of a logistic regression output in R. I used training data to make the model.   * How do I test the logistic regression model developed on the training data on the data left out? My naive guess is to create a function then run each test same through that (not even sure how do pull that) but I have to imagine there's a better way.


title: CODA gleman.diag, Error in chol.default(W):

text: I'd like to use gelman.diag for an MCMC chain I ran in JAGS. It is very large, so I can't provide it. The chain contains several MVN distributions, and I use a wishart prior on the precision matrix. Just wondering what potential sources of this error might be? Thanks you


title: Explanation for R-squared as ratio of covariances and variances

text: I have code that calculates $R^2$ with summations $$R^2 = \frac{(\sum xy - \frac1n \sum x \sum y)^2}{(\sum x^2 - \frac1n \sum x \sum x) (\sum y^2 - \frac1n \sum y \sum y)},$$ which is equivalent to $$R^2 = \frac{cov(x, y) \cdot cov(x, y)}{var(x) \cdot var(y)}.$$ I know the code is correct by benchmarking, but I have never seen this form. Can someone please explain or provide a reference? Thanks! FWIW, the code is built for speed. It does rolling regressions and can quickly find each summation by differencing a cumulative sum.


title: Managing the experiment loop - experiment organization tool

text: Here is loop of actions I do very often in my machine learning research:   1. Run an experiment with a certain configuration (a certain classifier with certain settings, etc.).   2. Look at the log file. Try to think what configuration to change in order to improve the results.   3. Rename the log file of the experiment, so that I remember what configuration it came from. For example, from "experiment.log" to "experiment.bayesian.log" or "experiment.decisiontree.log".   4. Change the configuration and return to step 1. This loop has many technical actions, and many opportunities to make a mistake. For example: I renamed the log, but didn't include one of the configuration params because I thought it is not significant. Then the results proved that it is significant, and I have to go over all the logs and remane them. For example, rename "experiment.bayesian.withoutboosting.log", and the new log file to "experiment.bayesian.withboosting.log" So, I wonder if there is a system that automates this process? Preferrably, I am looking for a generic tool, that does not depend on the programming language I use for my classifiers. Maybe a shell script, that changes a configuration file, runs a program, extracts the results from the log, and keeps them in a table. Do you know of anything like this?


title: Probability of data point being from distribution in normal mixtures

text: Using JMP, I was able to fit a distribution to a set of data, using the normal-2 mixtures model. It returns location (or mean), dispersion (standard deviation) and probability for each of the two normal distributions used to create the normal-2 mixtures. Now, I want to be able to take any data point from that population, and figure out the odds that data point came from each of the two distributions. Is there a way to do this?


title: SEM & Identification

text: Consider the following demand and supply model:   * Demand equation: $q=a_1p+a_2y+e^d$    * Supply equation: $q=b_1+p+e^s$  The supply equation is identified because $y$ can be used as an instrument for $p$. The demand equation is **not** identified as it stands. I am struggling to understand how **linear restrictions on parameters** can be used to identify the demand equation. For instance: 1) Set $a_2=0.5$ and use $y$ as an instrument in $(q-0.5y)=a_1p+e^d$ 2) Set $a_2=-0.5a_1$ and use $y$ as an instrument in $q=a_1(p-0.5y)+e^d$ **Why do 1) and 2) identify the demand equation?** Thank-you!


title: Methods for randomly allocating people between active and placebo in clinical trials

text: I'm looking for a way to allocate people to either group A or B. Data for the study group has 3 main variables:   1. Age 65-90   2. Male or Female   3. APOE + or APOE – (this is a blood test) I then need to allocate them to either A (active tablet) or B (placebo). Ideally, I need both groups A and B to be comparative, so fairly even. Total number of participants will be 120. Can anyone suggest a random allocation method for doing this? Or is it perhaps better to just randomly assign and completely ignore the above variables?


title: Grouping cases on certain criterion

text: You experts have been helpful so many times that I read others' threads, but this time I have to ask a new question. Long story short I need to _group all my cases into 4 categories_. Could someone please help me figuring out the procedure in spss? (I have individuals' scores on 2 variables - fashion innovativeness and fashion opinion leadership - and other data which is irrelevant to my question) _Fashion innovators_ \- persons scoring more than one standard deviation above the mean on fashion innovativeness but less than one SD above the mean on fashion opinion leadership. _Fashion opinion leaders_ \- persons scoring more than one SD above the mean on fashion opinion leadership but less than one SD above the mean on fashion innovativeness. _Innovative communicators_ \- persons scoring more than one SD above the mean on both fashion innovativeness and fashion opinion leadership. _Fashion followers_ \- persons scoring less than one SL above the mean on both fashion innovativeness and fashion opinion leadership. How do I create groups with these conditions? I need that in order to give a demographic profile of each group and compare them on certain traits. Sorry, I am truly on bad terms with spss. Thanks much for you help!


title: Predicting football match winner based only on the outcome of previous matches between the two teams

text: I'm a huge football (soccer) fan and interested in machine learning too. As a project for my ML course I'm trying to build a model that would predict the chance of winning for the home team, given the names of the home and away team. (I query my dataset and accordingly create datapoints based on previous matches between those 2 teams) I have data for several seasons for all teams however I have the following issues that I would like some advice with. The EPL (English Premier League) has 20 teams which play each other at home and away (380 total games in a season). Thus, each season, any two teams play each other only twice. I have data for the past 10+ years, resulting in 2*10=20 datapoints for the two teams. However I do not want to go past 3 years since I believe teams change quite considerably over time (ManCity, Liverpool) and this would only introduce more error into the system. So this results in just around 6-8 data points for each pair of team. However, I do have several features (up to 20+) for each data point like full-time goals, half time goals, passes, shots, yellows, reds, etc. for both teams so I can include features like recent form, recent home form, recent away form etc. However the idea of just having only 6-8 datapoints to train with seems incorrect to me. Any thoughts on how I could counter this problem? (if this is a problem in the first place)


title: Is there an efficient way to discriminate space based on K-Means results?

text: Suppose we done K-Means and got K centroids of clusters and we want to tag new points based on those K centroids. **UPDATE:** These K centroids are given to me, so I can't go for another clustering algorithm. Also I have a large number of 2d points and a big k, so it's very important to get appropriate complexity. The solution came into my mind is to do:               For every p in newPoints do:         For every c in centroids do:             calculate distance between c and p             if distance < minDistance                 minDistance = distance                 p.tag = c.tag         end for     end for      But the complexity of this solution is O(K*N) where N is the number of new points. I would like to know if there is a solution with less complexity.


title: hypothetical measure of variability similar to entropy

text: Variance has the following properties:   1. $Var(cX)=c^2Var(X)$   2. For independent variables $Var(X+Y)=Var(X)+Var(Y)$. Range of a rv has the following properties:   1. $Range(cX)=|c| Range(X)$   2. For independent variables $Range(X+Y)=Range(X)+Range(Y)$. Entropy has the following properties:   1. $H(cX)=H(X)$   2. For independent variables $H(X,Y)=H(X)+H(Y)$. My main question: Does there exist a _characteristic_ of a random variable such that:   1. $New(cX)=New(X)$   2. For independent variables $New(X+Y)=New(X)+New(Y)$. ? (the only difference is "plus" instead of a vector random variable) Maybe something can be said if the first property is replaced by $New(cX)=|c|^3New(X)$?


title: Degrees of freedom ANOVA table for regression

text: I've calculated the following ANOVA table for a linear regression in R:               anova(lm(mpg ~ drat, mtcars))                Analysis of Variance Table          Response: mpg      Df Sum Sq Mean Sq F value Pr(>F)      drat 1 522.48 522.48 25.97 1.776e-05 ***      Residuals 30 603.57 20.12      ---       I'm wondering how the degrees of freedom for the residuals are calculated. I can see it is n-2, but why minus 2?


title: Time series anomalies

text: Suppose I have some time series. I already removed seasonality, so it consists of trend, some events and error. I need to find the events, their length in periods,and their strength they had on the trend. My idea is to divide data into several "windows", and calculate the length and strength Do you know of a similar technique?


title: How to understand the variance in this setting?

text: All - I have a dataset originated from a practical setting, but it is not clear to me how to interpret it. Let me try to frame the setting and see if it makes sense: A manufacturer produces a type of product that use noise level (dB) as one of its quality metrics. All products during a production cycle are expected to have similar or consistent noise levels. Further, Let's say this "consistency requirement" is defined here as "the max and min noise level should not differ more than $x$ percentage point from the mean". My question are two-fold:   * Is it possible to use some distribution metrics, for example, standard devision, to make claim such as: if the $sd < a$, then it will meet the consistency requirement?   * It seems intuitively that the more products produced in the cycle, the harder it is to meet the consistency requirement. Statistically, each product's noise is independent variable, the larger the population, the larger the variation of sum - is this a correct understanding/statement? Thanks Oliver


title: What exactly does the 'boxcox' function in R do?

text: I am familiar with the power transform family and I know how to estimate the MLE for $\lambda$ for given samples of a random variable. I have been using the 'boxcox' function in R for a sample of a random variable that depends linearly on a latent variable (i.e. a linear model) and it works very well. The problem is that I don't quite know how it works on such linear models. The official reference won't give clear clues on how to answer my question. Can you provide some information on the answer to my question?


title: Interpreting the changing significance of a variable in a Logistic Regression

text: I am testing the relationship between two variables (Similarity index; Years patenting) on a binomial dependent variable via a Logistic regression in SPSS. I tried two models as you can see bellow in the links. In the first model I tested, only one variable (Years patenting) was shown significant. In the second model, where I added the interaction between the two variables, that interaction became significant as well as the other variable (Similarity index). The variable that was shown significant in the first model (Years patenting) was not significant anymore. How would you interpret these results? Model 1: ![enter image description here](http://i.stack.imgur.com/pRkG9.png) Model 2: ![enter image description here](http://i.stack.imgur.com/JnhB3.png) Also - how would I know the direction of the relationship between the predictor and the outcome? Do I look at the sign for the regression coefficient B, or do I look if Exp(B) is larger or smaller than 1?


title: Confidence Interval for Inverse Gamma Distribution

text: I would like to understand if there exists any method to find confidence interval for the parameters of inverse gamma distribution. Thanks in advance!


title: How to report results from eight paired t-tests in a scientific paper?

text: What would be the most efficient way to report test results from several paired t-tests in a scientific paper? If I gave a bar graph per test result (something like this) would it be too much? Should I do a table instead? Or both? Give graphs only for the 'best' results?


title: Cluster-randomized trial: mixed model hypothesis testing

text: I want to estimate the effect of randomly assigned intervention. The outcome is measured at the individual level, but the individuals are assigned to groups which influence eachother a lot, and it is the groups which are assigned to treatment or control. I need to test the null hypothesis that the intervention had no effect. I think this is a case where I want to estimate random effects at the group level and a fixed effect for the treatment (I've used lmer before), but I'm not totally sure, and even if that's right, I'm not sure how to take the next step from there to test the null I want to test.


title: Bootstrapping in SAS - PROC LOGISTIC - Next steps ? how to score / perform diagnostics?

text: My question is as follows. I am referencing the following paper by David Cassell - wherein David talks about bootstrapping techniques in SAS using PROC SURVEYSELECT (many thanks to David - truly a seminal piece of work): http://www2.sas.com/proceedings/forum2007/183-2007.pdf Essentially, David boostraps "n" replicates and uses a PROC MEANS to arrive at the average of the parameter of interest. I used the approach to perform PROC LOGISTIC --> that is, I performed 2000 regressions, averaged the parameter estimates and now I have the "final" model. But I need to score other datasets using this model. I also need to perform regression diagnostics using this (final) model. My questions are as follows:   1. How do I now take these final estimates (betas) and the modify the INMODEL= or the OUTMODEL= dataset produced/consumed by PROC LOGISTIC to score my new/training/any other dataset ? OR Can I "force" PROC SCORE to consume these bootstrap estimates and then - score another dataset ? [From my research on the internet, apparently modifying the INMODEL or OUTMODEL datasets is frowned upon !]   2. Secondly, how do I perform regression diagnostics - outliers, influential observations with these "new" bootstrap estimates. In some sense, I don't want to manually code the regression equation to get at my predicted values, compute the studentized residuals, the hat matrix (leverage), DFBETAS, Cook's D etc if you will.  While I am very happy that I bootrapped my regressions and proved to my team that, sampling variability does indeed play a crucial part in parameter estimation in small samples, I am only about 50% finished with the entire gamut of modeling. Any help would be greatly appreciated ? How do my fellow SAS bootstrap modelers go about SCORING/DIAGNOSING/REMEDYING their models once the bootstrap estimates have been estimated. Are David or Peter Flom listening :-) ??? Is Lex Jansen actually a real person or an internet moniker ??? :-).....Shout out to Rick Wicklin :-) ! Thanks Sincerely R


title: Area under the ROC curve or area under the PR curve for imbalanced data?

text: I have some doubts about which performance measure to use, area under the ROC curve (TPR as a function of FPR) or area under the precision-recall curve (precision as a function of recall). My data is imbalanced, i.e., the number of negative instances is much larger than positive instances. I am using the output prediction of weka, a sample is:               inst#,actual,predicted,prediction     1,2:0,2:0,0.873     2,2:0,2:0,0.972     3,2:0,2:0,0.97     4,2:0,2:0,0.97     5,2:0,2:0,0.97     6,2:0,2:0,0.896     7,2:0,2:0,0.973      And I am using pROC and ROCR r libraries.


title: SPSS Contrasts for a Five-Way Within-Subjects ANOVA (Yes, You read correctly)

text: I have a 4 x 4 x 2 x 2 x 2 repeated measures design and I'm trying to test a contrast for just the first variable (here called "SOA"), averaging over the others, in SPSS. Initially, I thought it would simply be entered into syntax as: /WSFACTOR=SOA 4 special(1 1 1 1, 1 -3 1 1) But I get back a message reading "Number of specified values for the SPECIAL contrast matrix is not the same as squared number of levels. Execution of this command stops." Just wondering what the problem might be? Thanks.


title: How can we combine learnings from multiple experiments in a single causal model?

text: I would like to use a causal network modelling to model the interaction of several variables and the effects of interventions. I have measurements for all priors of the model, that is without any interventions, and a well defined structure as DAG. Also, assume there are no hidden variables. We have the possibility to run experiments and measure the causal effects of interventions on specific variables. Assume we can set one variable per experiment. That is we have measurements of the type $P(A / do(B) )$ using Pearls’ notation. Different experiments lead to different and new learnings and measurements of the type $P(A / do(C) )$. Is there a framework on how to combine these learnings from many different experiments in order to infer quantities such as $P(A / do(B), do(C) )$ without running an experiment that sets both $B$ and $C$?


title: Why do we need $R^2$?

text: In linear regression, the $R^2$ value is the square of the correlation between predicted values and observed values. But why do we need the $R^2$ value? Why not just use the correlation coefficient? Just like the correlation coefficient, $R^2$ is scale-less (i.e. values are always between 0 and 1), so I can't see why there's a need for $R^2$. I would imagine it is something to do with the fact that the correlation coefficient can be negative, but don't really see why this would be a problem.


title: Testing the hypothesis on clustering

text: I have a number of samples. For each, there is a time course of multivariate data defined, with $t$ timepoints ($t < 50$) and $n$ variables ($n > 100$). We have noted that the time courses of a certain further variable $X$ _seem_ to split this group of samples in two. We would like to use the $n \times t$ variables to test this hypothesis. How would you handle this problem? I'm a bit lost at the moment. I was considering using PCA or MPCA and than an acceptable clustering technique to see whether we can recover the original clustering. However, several alternative clusterings could be possible in the given data set (for example, we can have red and blue circles and triangles). I'm not sure how to test the following statement: There exists a _substantial_ (significant? statistically significant?) set of variables that show the same pattern of clustering as $X$.


title: Logistic Regression with different priors

text: I am using standard logistic regression for classification with reasonable results. As expected I get a probability of 0.5 for query points "far away" from the data. However I would like to assign these points to one of the classes (or to a third class if that is easier). Is this possible?


title: Constrained fixed effects

text: I need to show that constraining the sum of group level fixed effects (in this case, zero) has no effect on the coefficients of the regressors. My intuition is that each constrained d_i is a perfect linear transformation of the non constrained constants. I looked in Greene, and Cameron and Trivedi, but I did not find a formal treatment of it.


title: Returning the inverse of a matrix in a quadratic program (SVM) in cvx optimization package

text: I am solving the dual QP of an SVM, and using the RBF kernel. As you know, the objective function is of the form $$f(\alpha) = \alpha^T Q \alpha $$ where $\alpha$ is the optimization variable and $Q$ is some positive semidefinite matrix. When $Q$ is dense, I wish to see what cvx (or any other optimization packages if you can point me to one) uses as $Q^{-1}$ when solving the QP. Any ideas?


title: plotting a main effect in the context of a significant interaction

text: I have a model that has two significant main effects and where the interaction between the two variables (age and %gray) is also significant in a multiple regression. When the interaction (age * %gray) is present in the model, the beta of age changes sign with respect to the model where the main effects alone are included. Is it possible to plot the data (independent variable vs. dependent) in a way that reflects the new sign of the main effect, after the interaction term has been included?


title: Normalize row or column while each row is an observation

text: Suppose I have a matrix compose of row as each observation, column as each property and I want to calculate the distance between each observation. In this case I think I should normalize each column, instead of row, to a unit vector? I think that would make more sense since each column can have different unit. In contrast if I normalize it by rows to unit vectors, it will be dominated by the property that has larger average value.


title: Using regression to prove that X is a necessary / condition condition for Y

text: I have been thinking about how to use regression to prove whether X factor is a necessary or a sufficient condition for Y. I am employing these terms in the traditionally logical sense, i.e.   * If no X, then no Y --> X is necessary for Y   * If X, then Y --> X is sufficient for Y The reason I'm asking is because in social sciences theories we talk about necessary and sufficient conditions all the times (e.g., "(If) No bourgeoisie, (then) no democracy"). At the same time, the most prevalent tool is regression, yet it seems ill-fitted to investigate necessary / sufficient causes. Indeed, rarely do we talk about necessity and sufficiency in the context of regression. The usual interpretation is that X "causes" Y in an additive sense -- i.e. a certain change in X is associated / caused a certain change in Y. So, is it possible to interpret regression results as evidence of necessity versus sufficiency? In contrast, a simple 2-by-2 table of X and Y would make the necessity / sufficiency relationship very clear. P/S: the motivation needs not be restricted to social sciences. For example, if I have a dataset with DV as "Fire occurring", and IVs as presence of fuel, presence of match, presence of lightning, etc. Can we use regression to prove that fuel is necessary but insufficient for fire?


title: recursive feature elemination in R with caret

text: i work with R `caret` software package to select the most important features from some set of data. My response is a factor of multiple classes (e.g. nominal Attribute) and my features are nominal as well as numerical. I read in the software manual that `rfe()` works as following > For classification, ROC curve analysis is conducted on each predictor. For > two class problems, a series of cutoffs is applied to the predictor data to > predict the class. The sensitivity and specificity are computed for each > cutoff and the ROC curve is computed. The trapezoidal rule is used to > compute the area under the ROC curve. This area is used as the measure of > variable importance. For multi-class outcomes, the problem is decomposed > into all pair-wise problems and the area under the curve is calculated for > each class pair (i.e. class 1 vs. class 2, class 2 vs. class 3 etc.). For a > specific class, the maximum area under the curve across the relevant pair- > wise AUC's is used as the variable importance measure. So far so good. But i do not understand how the ROC curve is calculated if my feature is nominal. For illustration i have some toy data set up and perform `rfe()`. How is the AUC for `Geschlecht` calculated??               rm(list=ls())          library(klaR)     library(caret)     library(matlab)     library(FSelector)          Faktor<-factor((sample( LETTERS[1:4], 10000, replace=TRUE, prob=c(0.1, 0.2, 0.65,   0.05))))     alter<-abs(rnorm(10000,30,5))     HF<-abs(rnorm(10000,1000,200))     Diffalq<-rnorm(10000)     Geschlecht<-sample(c("Mann","Frau", "Firma"),10000,replace=TRUE)     data<-data.frame(Faktor,alter,HF,Diffalq,Geschlecht)               set.seed(5678)     flds<-createFolds(data$Faktor, 10)          train<-data[-flds$Fold01 ,]         test<-data[flds$Fold01 ,]               features <- c("HF","alter","Diffalq", "Geschlecht")          formel<-as.formula(paste("Faktor ~ ", paste(features, collapse= "+")))          nb<-NaiveBayes(formel, train, usekernel=TRUE)          pred<-predict(nb,test)          test$Prognose<-as.factor(pred$class)          start<-proc.time()     set.seed(12345)     nbProfile<-rfe(Faktor ~ alter + HF + Diffalq + Geschlecht, train, sizes=1:4,repeats=3,  rfeControl =  rfeControl(functions = nbFuncs, method = "cv"))     ende<-proc.time()     zeit<-ende-start          nbProfile      Thanks in advance!


title: Practical Implementation of Speech Recognition HMM

text: I'm trying to implement a GMM/HMM for phoneme recognition where I have for each phoneme a 3-state left to right HMM model with start and end states with no emission. The emission probabilities are modeled with Gaussian mixtures. I'd like to build an HMM architecture as proposed in paper "Speaker- independent phone recognition using hidden Markov models" by Kai-Fu Lee Figure 3-3 (a). That is to have a start state and end state with no emission and in between to have all the phoneme hmms in parallel for continuous phoneme recognition. Is there any reference or a paper that indicates how to practically implement this or how to implement the start and end states without emitting symbol?   When I refer to implement, I mean to implement Baum-Welch, Viterbi, etc. Thanks!


title: Formula for non-linear regression in R

text: I want to know if it is possible for a library in R to evaluate the association of independent variables and create a formula? I am trying to come up with a model to predict power consumption of a machine, using some hardware counters and performance attributes. When I use linear regression, I have no problem since I could represent my formula like `power~lm(a1+a2+a3+a4)`, but for the non-linear case, I am not sure what would be the formula or which model should I choose. I would want to have a way to do this:               power ~ <some-non-linear-reg-pkg>(a1+a2+a3+non-linear(a4))      I reviewed some packages for non-linear regression such as `nls` and `gnm`, and they expect a formula to be provided by the user. I am however able to identify which variables have linear associations and which are non-linear (by performing correlation tests), the problem is building a formula out of them.


title: On the connection between SSE and absolute deviation from the centroids

text: Is there any connection between sum of squared error SSE and the absolute deviation from the centroids after clustering. More formally, I have clustered $T=\\{x_i\\}, i\in\\{1,\ldots,n\\}$ and the results are $c$ clusters: $T^C=\\{C_j\\}, j\in\\{1,\ldots,c\\}$ (As a side note, I appreciate your comments on the notations). All records $x_i$ are assigned to clusters $G_j$ with the centroids $C_j=MEAN(x_i), x_i\in G_j$, so SSE is calculated as: $\sum\limits_{\underset{x_i\in G_j}{i=1}}^{n}(x_i-C_j)^2$ I am seeking a way to calculate/estimate the following: $\sum\limits_{\underset{x_i\in G_j}{i=1}}^{n}|x_i-C_j|$ However, I am not sure whether there is such a connection or not. My algorithm is something similar to K-means and Euclidean distances are used. Thanks


title: Consequences of modeling a non-stationary process using ARMA?

text: I understand we should use ARIMA for modelling a non-stationary time series. Also, everything I read says ARMA should only be used for stationary time series. What I'm trying to understand is, what happens in practice when misclassifying a model, and assuming `d = 0` for a time series that's non-stationary? For example:               controlData <- arima.sim(list(order = c(1,1,1), ar = .5, ma = .5), n = 44)      control data looks like this:                [1]   0.0000000   0.1240838  -1.4544087  -3.1943094  -5.6205257      [6]  -8.5636126 -10.1573548  -9.2822666 -10.0174493 -11.0105225     [11] -11.4726127 -13.8827001 -16.6040541 -19.1966633 -22.0543414     [16] -24.8542959 -25.2883155 -23.6519271 -21.8270981 -21.4351267     [21] -22.6155812 -21.9189036 -20.2064343 -18.2516852 -15.5822178     [26] -13.2248230 -13.4220158 -13.8823855 -14.6122867 -16.4143756     [31] -16.8726071 -15.8499558 -14.0805114 -11.4016515  -9.3330560     [36]  -7.5676563  -6.3691600  -6.8471371  -7.5982880  -8.9692152     [41] -10.6733419 -11.6865440 -12.2503202 -13.5314306 -13.4654890      Assuming I didn't know the data was `ARIMA(1,1,1)`, I might have a look at `pacf(controlData)`. ![pacf\(controlData\)](http://i.stack.imgur.com/IOXJf.jpg) Then I use Dickey-Fuller to see if the data is non-stationary:               require('tseries')     adf.test(controlData)          # Augmented Dickey-Fuller Test     #     # data:  controlData     # Dickey-Fuller = -2.4133, Lag order = 3, p-value = 0.4099     # alternative hypothesis: stationary          adf.test(controlData, k = 1)          # Augmented Dickey-Fuller Test     #     #data:  controlData     # Dickey-Fuller = -3.1469, Lag order = 1, p-value = 0.1188     # alternative hypothesis: stationary      So, I might assume the data is ARIMA(2,0,*) Then use `auto.arima(controlData)` to try to get a best fit?               require('forecast')     naiveFit <- auto.arima(controlData)     navifeFit     # Series: controlData      # ARIMA(2,0,1) with non-zero mean      #      # Coefficients:     #          ar1      ar2     ma1  intercept     #      1.4985  -0.5637  0.6427   -11.8690     # s.e.  0.1508   0.1546  0.1912     3.2647     #     # sigma^2 estimated as 0.8936:  log likelihood=-64.01     # AIC=138.02   AICc=139.56   BIC=147.05      So, even though the past and future data is ARIMA(1,1,1), I might be tempted to classify it as ARIMA(2,0,1). `tsdata(auto.arima(controlData))` looks good too. Here is what an informed modeler would find:               informedFit <- arima(controlData, order = c(1,1,1))     # informedFit     # Series: controlData      # ARIMA(1,1,1)                         #     # Coefficients:     #          ar1     ma1     #       0.4936  0.6859     # s.e.  0.1564  0.1764     #     # sigma^2 estimated as 0.9571:  log likelihood=-62.22     # AIC=130.44   AICc=131.04   BIC=135.79      1) Why are these information criterion better than the model selected by `auto.arima(controlData)`? Now, I just graphically compare the real data, and the 2 models:               plot(controlData)     lines(fitted(naiveFit), col = "red")     lines(fitted(informedFit), col = "blue")      ![tsPlots](http://i.stack.imgur.com/sy3YR.jpg) 2) Playing devil's advocate, what kind of consequences would I pay by using an ARIMA(2, 0, 1) as a model? What are the risks of this error? 3) I'm mostly concerned about any implications for multi-period forward predictions. I assume they would be less accurate? I'm just looking for some proof. 4) Would you suggest an alternative method for model selection? Are there any problems with my reasoning as an "uninformed" modeler? I'm really curious what are the other consequences of this kind of misclassification. I've been looking for some sources and just couldn't find anything. All the literature I could find only touches on this subject, instead just stating the data should be stationary before performing ARMA, and if it's non-stationary, then it needs to be differenced d times. Thanks!


title: Difference between Excel's exponential chart trend line and GROWTH/LINEST function?

text: I have a exponential trendline that was generated by excel ("Add Trendline") on a chart. I've plugged in the equation (from the "Display Equation on Chart" option, up to 14 decimals) into a spreadsheet and compared the projected values with those values projected by the GROWTH and LINEST functions. In my comparison, I noticed that the interpolated values differ. Can anyone explain how the methods excel uses to derive these trendline equations differ? Is one method more appropriate than the other? Thank you.


title: Can AUC decrease with additional variables?

text: I'm fitting a logistic regression model to predict probabilities from a set of variables. I'm comparing two such models, say `M1` and `M2`. The only difference is that `M2` includes all the variables of `M1` plus a few more variables. The idea is to see which variables are useful in predicting my dependent variable. I expected that AUCs should be non-decreasing with the addition of new variables. If the new variables have predictive power, they should increase the AUC, if they don't, then the AUC should be unaffected. But I find that AUC actually decreases as I add a particular set of new variables. What could be the issue here? I'm using `predict()` to get the predicted probabilities. Does it automatically drop all the statistically insignificant variables when calculating the predicted value? Could this be the cause of the drop in AUC?


title: How to down weight correlations in my microarray analysis?

text: **Background:** I have been tasked in one part of my analysis to reproduce a method used in another study as follows in bullet points form:   * Microarray data from a number of time points   * Calculate correlation between two subsets of genes, set 1 and set 2   * On the basis of a KS test they determine there is a signal in the data suggesting that set 1 and set 2 are higher correlated than by chance   * To down weight spurious correlations, they randomize all the time points 10^5 times and calculate r between each pair of genes each time - taking "p values" as the number of times a pair of genes with randomized time points is correlated at least as well as the real pair with ordered time points (divided by 10^5). My Problem: I have fewer timepoints, with each measurement having only 24 permutations of its timepoints possible. Hence, step 4 in which measurements have their time points randomized is problematic to reproduce in any way that is useful. Hence I am considering predicting new time points which will be used in calculating correlations, thereby allowing me to carry out their "empirical" approach to assigning p-values. (Note: they dont really treat them as p-values, just scores for down-weighing spurious correlations). **Question:** What methods are there of down-weighting spurious correlations? On the basis of the above filler, would predicting new time points for my data be a complete statistical disaster? I am sure its not the best idea but the exact reasons for why that is are eluding me.


title: What's a good way to mentally visualize n dimensions in a k means

text: I've been using k-means to do some clustering and one of the ideas I'm struggling with is the n dimensions aspect. If I were clustering housing prices vs sq. feet its just a simple 2d graph. That I can visualize and it makes sense. Where it gets difficult is when it's more than 3 features. If I understand properly, each feature is plotted in it's own dimension. How does something with 3+n dimensions look, conceptually? What's a good way for me to think about how the data is represented? I know this question is a bit vague (which is a result of my lack of understanding in this area), but any help is appreciated.


title: Ftting a mixture of two Gaussians

text: I want to fit a mixture of two gaussian densities to my financial data. The data can be found here: http://uploadeasy.net/upload/2a7mw.rar the variable is called dat. The probability density of a mixture is given by: \begin{align} f(l)=\pi \phi(l;\mu_1,\sigma^2_1)+(1-\pi)\phi(l;\mu_2,\sigma^2_2) \end{align} The quantile can be computed by using a numerical algorithm to solve the following: \begin{align} \alpha=P(L \leq VaR_\alpha) = \pi F_1(Quantile_\alpha;\mu_1,\sigma^2_1)+(1-\pi) F_2(Quantile_\alpha;\mu_2,\sigma^2_2) \end{align} I use mixtools in R:               install.packages("mixtools")     library(mixtools)     mixture<-normalmixEM(dat,k=2,fast=TRUE)      This uses the EM algorithm. I now want to calculate the 0,95 quantile of the mixture distribution. I do a loop, a kind of a grid search, I assume, that the quantile (due to the characteristics of my data) will be below 0.3. So the loop ends at 0.3               pi<-mixture$lambda[1]         mu1<-mixture$mu[1]     mu2<-mixture$mu[2]         sigma1<-mixture$sigma[1]     sigma2<-mixture$sigma[2]          quantile<-0     probabilitylevel<-0.95     dummy1<-0          # the loop lasts for about 20-40 seconds     for (i in 1:100000){     quantile[i]<-i/(1000000/3)     }     dummy1<- probabilitylevel - ( pi * pnorm(quantile,mean=mu1,sd=sigma1) + (1-pi) * pnorm(quantile,mean=mu2,sd=sigma2))          min(abs(dummy1))     which.min(abs(dummy1))     quantileresult<-which.min(abs(dummy1))/(1000000/3)      the result               quantileresult      is `0.025371` which seems to be correct, if control it with:               pi * pnorm(quantileresult,mean=mu1,sd=sigma1) + (1-pi) * pnorm(quantileresult,mean=mu2,sd=sigma2)      I look at the plot:               plot(density(dat),col="red")     curve(expr=pi*dnorm(x,mu1,sigma1)+(1-pi)*dnorm(x,mu2,sigma2),lwd=2,col="black",add=TRUE)     curve(dnorm(x,mean(dat),sd(dat)),add=TRUE,lty=3,col="orange",lwd=2)      which gives ![quantilecorrect](http://i.stack.imgur.com/DAoHg.png) It looks like, that the mixture normal (black) is fitting the data way better. The dashed orange line is the univariate normal distribution fitted to the data set. It is fitting the data not as good as the mixture density, is this correct interpreted? Finally, we look at the single densities and compare it to the mixture:               plot(density(dat),col="red")     curve(dnorm(x,mu1,sigma1),add=TRUE,lty=2,col="darkgreen")     curve(dnorm(x,mu2,sigma2),add=TRUE,lty=2,col="blue")     curve(expr=pi*dnorm(x,mu1,sigma1)+(1-pi)*dnorm(x,mu2,sigma2),lwd=2,col="black",add=TRUE)      which gives the following plot: ![densitiescompared](http://i.stack.imgur.com/6dyKD.png) The first density has a higher peak, the second density is shifted to the left and has a lower peak, higher variance. Are my calculations and interpretations correct?


title: Uniform distribution - a simple question

text: I know this is a really simple one, but for some reason I see various ways and I'm not sure which one should I follow. So - I have $Y_1,...Y_n \sim U(1,3)$ and I want to know $P(y<c)$. The answer that I have for this (and it's different from how I see in Wiki) is: $c+\int_1^c 0.5\,dy $ Is this the correct way? If so, what's the logic behind it?


title: Calculating probability of discovery

text: I have a planner that can evaluate N arbitrary states, and calculate their fitness. The domains it evaluates have no explicit "end state", so it has an infinite horizon. What are good methods for calculating when it should stop evaluating states? My current approach is to estimate the probability of the next evaluation being the next "best" state, and to stop if this probability falls below a certain threshold. e.g. After evaluating each state, I increment a counter, and after a new best fitness if found, I reset that counter, resulting in a sequence like:               steps = [0, 1, 2, 0, 1, 2, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4]      I can then iterate over this and form a simple discrete probability of each step preceding the next best fitness:               from collections import defaultdict     counts = defaultdict(int)     totals = defaultdict(int)     for step,next in zip(steps, steps[1:]):         counts[step] += next == 0         totals[step] += 1     for k in sorted(counts.keys()):         print k,counts[k],totals[k],'%.2f' % (counts[k]/float(totals[k]),)      resulting in:               0 0.43     1 0.00     2 0.50     3 0.00     4 0.00     5 0.00     6 0.00     7 0.00     8 0.00     9 0.00     10 1.00      How do I take these discrete probabilities and form a general approximation for the probability of `P(next step is the best | current step since last best)` given an arbitrary step integer (especially for high step counts I may not have reached)?


title: Generating a 3-month half life weighting series in R

text: I am using the lm() and princomp() functions in R to perform regressions on foreign exchange time series. I would like to weight the regressions (and PCA) such that 50% of the influence on the regression comes from the past 3 months, 25% from the previous 3 months, etc, but in a smooth fashion. Both functions take such a weighting series. What is a simple formula for generating a decaying weights series with a half life of 3 months (or any other period, for that matter), in R? Ideally the sum of the weights will equal 1, although I don't think this is mandatory in lm().


title: Distribution of an observation-level Mahalanobis distance

text: If I have a multivariate normal i.i.d. sample $X_1, \ldots, X_n \sim N_p(\mu,\Sigma)$, and define $$d_i^2(b,A) = (X_i - b)' A^{-1} (X_i - b)$$ (which is sort of a Mahalanobis distance [squared] from a sample point to the vector $a$ using the matrix $A$ for weighting), **what is the distribution of** $d_i^2(\bar X,S)$ (Mahalanobis distance to the sample mean $\bar X$ using the sample covariance matrix $S$)? I am looking at a paper that claims it is $\chi^2_p$, but this is obviously wrong: the $\chi^2_p$ distribution would have been obtained for $d_i^2(\mu,\Sigma)$ using the (unknown) population mean vector and covariance matrix. When the sample analogues are plugged in, one ought to get a Hotelling $T^{\ 2}$ distribution, or a scaled $F(\cdot)$ distribution, or something like that, but not the $\chi^2_p$. I could not find the exact result either in Muirhead (2005), nor in Anderson (2003), nor in Mardia, Kent and Bibby (1979, 2003) . Apparently, these guys did not bother with outlier diagnostics, as the multivariate normal distribution is perfect and is easily obtained every time one collects multivariate data :-/. Things may be more complicated than that. The Hotelling $T^{\ 2}$ distribution result is based on assuming independence between the vector part and the matrix part; such independence holds for $\bar X$ and $S$, but it no longer holds for $X_i$ and $S$.


title: Compare coefficients from two separate panel regressions in Stata

text: I am trying to compare the coefficients of two panel data regressions with the same dependent variable. What I am aiming at is the following:                   y1 = c + β x         y2 = c + β x      In Stata                   xtreg y1 x i.z         xtreg y2 x i.z      I want to check whether the βs are significantly different. With two regular regressions I would use something like the following code in Stata to test a cross-equation restriction:                   sureg (y1 x ) (y2 x )         lincom [y1]x - [y2]x      However Stata is explaining that this is not possible when I try to use `xtreg`. I have searched a lot of different sites. However I do not have a clue, it only made me more confused. Hopefully you guys can help. ## edit I am using a Fixed effect models


title: Using principal components in a linear discriminant analysis for a diagnostic test

text: I am interested in building a linear discriminant function to discriminate between 2 groups, out of 60 variables. (I'm planning to select the most discriminative of the variables for a future diagnostic test.) I have calculated the area under the ROC curve for each of these variables individually and none have an AUC greater than 0.73. I have a fairly small sample of 50 healthy and 50 diseased individuals (these are the two groups). I have tried to reduce the number of variables using principal component analysis. There are 3 components accounting for 83% of the variation. But unfortunately, all 60 variables have similar weightings (loadings) in the 3 components, so I can't pick just few. I would ordinarily pick the highest weighted variables and then incorporate them in a linear discriminant function, but 60 is too many, especially given the small sample. I wondered if, **rather than** use the 60 variables, it is possible to use the 3 **principal components themselves** to build a linear discriminant function? I am using Stata for the analysis. Does anyone have any suggestions?


title: Estimating two-stage selection model with probit reduced form equation and ols structural form equation in R

text: I am trying to estimate a selection model of the form: $Z_i = 1[\alpha_0 + \alpha_1X_{1,i} + \alpha_2X_{2,i} + \delta_i$ > 0] $Y_i = \beta_0 + \beta_1X_{1,i} + Z_i + \epsilon_i$ where $1[]$ denotes the indicator function. The purpose of the model is to calculate the indirect effect of $X_1$ on $Y$ through $Z$, as well as the the direct effect. My first question is how to go about estimating this type of model, and how this estimation can be achieved in R. As far as I see it I have a few possible approaches: (1) Use a standard Heckman selection model, using OLS for both the reduced form and structural equations, using ivreg() in R. This will obviously ignore the constraint that $Z$ is bounded between 0 and 1. (2) Estimate the first stage with a probit model (i.e. $\delta_i \sim N(0,1)$), and the second stage using standard OLS. I understand that I could do this via manual 2SLS, but as far as I am aware the standard errors will be incorrect? Am I right in that this model is feasible, and if so, can you direct me to a method of achieving this in R? (3) Build a switching regression model (tobit-5) using the selection() function from sampleSelection package in R. I believe this model will estimate two equations for $Y$, one for where $Z_i=0$ and one where $Z_i=1$, and with a unique intercept and coefficients for each of the regressors in the outcome equations. The question then is how to get an estimate of the indirect effect of $X_1$ for each of these methods.   * If I use (1) or (2) then I imagine it might be possible to calculate the average marginal effect of $Z$ on $Y$, and the average marginal effect of $X_1$ on $Z$, then approximate the indirect effect by multiplying the two values?   * If (3) then could I take the fitted value under the estimated model for $Y$ where $Z=0$, and compare the mean to the mean of the fitted values under the estimated model for $Y$ where $Z=1$? This would then give me an estimate of the marginal effect of $Z$? Then use the same method as above and multiple this effect by the marginal effect of $X_1$ on $Z$? Many thanks in advance!


title: Nonlinear problem feature selection

text: I am trying to make prediction of students course score with neural networks. I have lots of parameters may affect to one course score like sex, GPA and many previous course scores. Before trying to make a neural network model I should make a feature selection, I guess. For example ABC440 may be affected ABC101, ABC 201 and ABC205 but I have 30 courses. Here I need to make a selection from other early courses. For a nonlinear problem like this, how can I make an efficient parameter selection? Can I use corelation analysis or linear regression analysis? Is it a correct way to use a linear analyse method for a nonlinear problem?


title: How can I avoid artificially inflating percentages?

text: I'm currently analyzing my data for one of my studies for my PhD dissertation in psychology. I'm hoping someone can help me with a problem I'm having... I've conducted a study in which participants play poker for one hour. I'm analyzing their decisions and giving every decision they made a code. If a decision had neutral or no expected value, the decision is coded as "2". If a decision has a negative expected value, I give a code of "0" for that decision. Lastly, if a decision has a positive expected value, I give it a code of "1". Participants vary in the number of hands they choose to play in the hour session. Some may play only 10 hands, while others may play 40 hands. Therefore, some individuals may have a large number of decisions coded, while others have much less. Here is my problem: I'm trying to determine how many, or what proportion of the participant's decisions have a negative expected value (i.e. poor decisions). So what I did was take the number of negative expected value decisions and divided it by total decisions. This creates a percentage of negative expected value decisions. The problem is that some people played very little hands and therefore had much fewer decisions. Having more poker experience is associated with playing fewer poker hands. Thus, for example, a more experienced individual may have made only 20 decisions in their playing session, while a less experienced playing may have made 50 decisions (due to having played more hands). If both players made 3 negative expected value decisions, then the percentage of negative expected value decisions for the player that played fewer hands would be 15% and the player that played more hands would have a percentage of 6%. Thus it appears as though individual with more experience are making poorer decisions! Playing fewer hands, and thus making fewer decisions, is artificially inflating the scores. Is there any way that I can remedy this? Perhaps by weighting or something of that nature? I work with SPSS. Any help would be greatly appreciated!


title: Is it possible to calculate the probability that two samples are drawn from the same distribution?

text: Given $X_1$..$X_n$ and $Y_1$..$Y_n$ drawn from unknown distributions $F(x)$ and $G(x)$ respectively, statistical tests such as two-sample Kolomogorov- Smirnov, Cramer-von Mises, and Anderson-Darling tests have been devised to test the null hypothesis $\mathcal{H}: F(x) = G(x)$ using various test statistics. But rather than hypothesis testing, I am more interested in quantifying the probability, $\mathcal{P}(X_1..Xn,Y_1..Y_n | F(x) = G(x))$ without knowing $F(x)$ or $G(x)$. Is this possible? The above tests give $\mathcal{P}(Z\ge{}z | F(x) = G(x))$ where $Z$ is the test statistic, but this is not what I want. I tried to find how the two-sample Kolomogorov test statistic was derived, but nothing useful came up... All the papers I have found either simply state the proposed test statistic or study its distribution under the null hypothesis. Given order statistics $X_{(1)}$..$X_{(n)}$ and $Y_{(1)}$..$Y_{(n)}$, is there anything useful we can say about $X_{(i)} - Y_{(i)}$ irrespective of the distributions $F(x)$ and $G(x)$?