title: Failing at linear regression / prediction on a real data set

text: I have a data set on which I'm trying to do regression, and failing. The situation:   * Thousands of battle robot operators are fighting battles among each other using battle robots.   * Some battle robots are strong and powerful, and others are weak; the strong ones win more often and deal more damage.   * Robot operators vary in skill, with the more skilled ones winning more often, and delivering more damage   * We have some summary information about the outcomes of their battles, but not all of the details.   * We know what battle robots they used in their battles, and how many times (including how many of those battles they won), and we know the total damage they dealt (of two kinds, damageA and damageB) in total   * Some robots are better at inflicting damageA, while others damageB   * For unknown battle robot operators based only on what robots they have used in battles (and how many times), we would like to estimate how much damage of each kind they would achieve, and what % of battles they have most likely won For example:   * John has used Robot A for 4 battles, and Robot B for 2 battles, and has dealt 240 units worth of DamageA   * James has used Robot A for 1 battle, and Robot B for 10 battles, and has dealt 1010 units worth of DamageA I can therefore estimate that Robot A probably deals 10 units of Damage A per battle, while Robot B deals 100 units of Damage A per battle, and thus if asked to estimate Damage A dealt by Matthew who has only played each of the two robots for 2 battles each, will estimate at 220 == (10*2 + 100*2). Unfortunately, the real data are not as clean and straightforward, probably because:   * There is a significant variance due to robot operator skill, e.g., a good operator could deal 20 units of damage with Robot A while a bad one only 5 units.    * There is some random variance due to opponents drawn in case of a small sample (e.g. somebody draws a strong opponent and loses despite having a better robot than the opponent), although eventually it would even out   * There may be some minor selection bias in that the best robot operators manage to pick the best robots to take into battle more often The real data set is available here (630k entries of known battle operator results): http://goo.gl/YAJp4O The data set is organized as follows, with one robot operator entry per row:   * Column 1 with no label - operator ID   * battles - total battles this operator has participated in   * victories - total battles this operator has won   * defeats - total battles this operator has lost   * damageA - total Damage A points inflicted   * damageB - total Damage B points inflicted   * 130 pairs of columns as follows:      * battles_[robotID] - battles using robot [robotID]     * victories_[robotID] - victories attained using robot [robotID] What I've done so far:   * Tried a couple of linear models using R `biglm` package which build a formula such as `damageA ~ 0 + battles_1501 + battles_4201 + ...` to try to get fitted "expected" values for each of the robots.   * Same, but removing the forced origin intercept by not including `0 +` in the formula   * Same, but also included the `victories_[robotID]` in the independent variables   * Same as before, but only selecting those robot operators whose victory numbers are close to their defeat numbers   * A linear regression model for `damageA ~ 0 + battles_1501 + battles_non_1501` where `battles_non_1501` are all the battles in robots other than robot model 1501. Then repeated for all the other robot types. I did sanity checks by looking at the predicted damageA and damageB values, as well as comparing the victories/battles ratio with the actual victories/battles ratio that we can actually precisely calculate for each of the robots. In all cases while the results weren't completely off, they were sufficiently off to see that the model isn't quite working. For example, some robots achieved negative damage numbers which shouldn't really happen as you cannot do negative damage in a battle. In case where I also used the known `victories_[robotID]` values in the formula, many of the `battle_[robotID]` coefficients ended up being somewhat large negative numbers, so I tried estimating for the "average" operator by `battle_[robotID] + victories_[robotID] / 2` but that also didn't give reasonable results. I'm somewhat out of ideas now.


title: Decision trees: nominal values unseen at training time

text: What is the best approach to deal with unexpected nominal values, unseen by decision trees at training time?


title: How Should I Find Statistically Significant Differences Between Two Markov Models?

text: Suppose I have N Markov Models of M states representing the behaviour patterns of 2 different groups (note: fully observable models, no hidden states), and have stored each model as a matrix of transition probabilities for a particular test subject. I would like to test whether the transition probability from any one state to any other is different in one group than the other. My first thought was to treat each of the MxM transition probabilities as a random variable, and test each hypothesis (i.e. that the means of the RVs differed) independently. The variable is a probability akin to one of the parameters of a multinomial distribution, so a binomial proportion-based confidence interval seems applicable. A correction for testing MxM hypothesises simultaneously from the same data also seems necessary. There are three questions that I'm left with now:   1. Is this the right approach at all?   2. Should I be using a different test?   3. Should I be applying a correction to the CIs for multiple hypothesis testing as described? I've read over this question already, but it's concerned with whether the matrices as a whole are identical, not whether any individual parameters differ. There doesn't seem to be much else available on SE for this topic.


title: If gauge charts are bad, why do cars have gauges?

text: It seems like data visualization experts generally disapprove of gauge charts (see here: What do you call a chart that looks like a half pie chart with a needle indicating a percentage?). The primary reason is that a gauge chart has a low data-to-ink ratio. Ever since I was exposed to these concepts (a few Tufte books), I generally agreed with them, but today it made me wonder: if gauges are so inefficient at communicating information, then why do cars/boats/planes have lots of gauges on their dashboards? And does the answer to _that_ question have some kind of relevance for creating software dashboards for large enterprises? Edited to include some additional information I found: I found a term, "glass cockpit", that refers to an airplance cockpit that has its mechanical gauges replaced with LCD screens. This gives credibility to the "convention" argument put forth by Wayne. http://en.wikipedia.org/wiki/Glass_cockpit Here's an iPad app that gives a dashboard-like readout of your car's telemetry, with no gauges to be seen. http://itunes.apple.com/us/app/dashcommand-obd-ii-gauge- dashboards/id321293183?mt=8 I also found a gross example of digital gauges for cars (viewer discretion is advised). http://www.chetcodigital.com/index-Automotive.htm


title: How to predict how long a worker will take to perform a task when performance on the task is only known for other workers?

text: I am trying to model a queue system (via simulation) to see if I can make different work assignments to improve the overall performance (based on a number of metrics). I have historical records of past tasks which includes when they entered the system, who was assigned to perform to the task, and the service time. In my simulator I want to have the option of replaying back these historical records to quantify how different assignment policies could have effected the performance measures in these counterfactuals. The workers perform at different speeds (i.e. have different service time distributions), so I would like to take this into account in my counterfactual simulations.   * **What is the best way to predict how longer worker $A$ will take on task $X$ when all I know about task $X$ is how long it took worker $B$ to complete it?**   * In other words, I have two distributions of the same type but with different parameters and given $x$ came from distribution $A$, **what would $x'$ be if it instead came from distribution $B$?** My current approach is to find the probability of $x$ in $A$, and then get a number for $x'$ that has similar probability in $B$. As an aside, I am very new to simulators and queuing theory in general so feel free to provide additional advice or links to resources if you think I could approach this a better way.


title: Neural network training on unlimited theoretical data

text: I'm considering using a neural network on financial time series but rather than train the network on actual data I am going to train on a model of the data which is perturbed by random noise. This being the case, I will potentially have unlimited amounts of model training data. However I don't want to generate a huge amount of data and then train the model as it might take a very long time to reach a solution, and I have no idea how much data to actually generate. What I am thinking of doing is training on a small amount of model data (say 5000 examples) and get the values for the hidden nodes, record these values, and then repeat again, thereby building up a distribution of values for each node. These distributions could then be bootstrapped to get the mean value per node, and the whole process would stop once the change in the bootstrapped mean values falls below a given threshold. Edit - the purpose of the network will be to classify/label the time series over the recent past as being in one of a finite number of states, e.g. trending up/down, moving sideways, in congestion etc. These states can be modelled using synthetic data with known labels, and then on real data the network's job will be to identify which state the real data most closely resembles. This will be used as input to a separate decision making process. My question is - is there any reason why this _would not_ be a valid approach to take?


title: Model without endogeneity correction has lower AIC than one with correction

text: I have two models, one with endogeneity correction (includes correction terms obtained using Heckman) and one without. The correction terms are significant in the second stage model, yet the AIC/BIC values of the model without the endogeneity correction is lower. Significant correction terms suggest presence of endogeneity and hence the need to correct for the same. However, the AIC/BIC values suggest that correction is not required. The model with the correction term makes more sense to me, but I would like to hear from the community about this anomaly and how do I justify using the model with the correction terms despite a higher AIC.


title: Number of independent samples for weighted samples?

text: I am taking random samples from one distribution, $f(x)$, but trying to get information about another distribution, $g(x)$. I have a weighting function, $w(x)=Cg(x)/f(x)$, to correct for this. The result is that I have $N$ independent samples, with different weights, $w_i$, attached to them. The question is then, what is a good estimate for the number of independent samples I really have. For example, if my weights are {0.49, 0.48, 0.01, 0.01, 0.01} then I have pretty close to 2 independent samples. If they are {0.3, 0.3, 0.4} then I have about 3. Presumably there is a quantitative way to do this. Also, how could I determine, given $f(x)$ and $w(x)$, what the efficiency of sampling is (i.e. How many independent samples of $g(x)$ do I, on average, get for $N$ samples of $f(x)$)?


title: Hilbert spaces and time series

text: Suppose that $\\{Y_{t}: t \in \mathbb{Z} \\}$ is a stationary zero mean time series. Consider the Hilbert space $\mathcal{H}$ generated by the random variables $\\{Y_t: t \in \mathbb{Z} \\}$ with inner product $$ \langle X, Y \rangle = E(XY)$$ and norm $$||X||^2 = E|X|^2$$ Consider the subspace $\mathcal{M}$ generated by the random variables $\\{Y_u: u \leq t \\}$. Why are future values found by projecting onto the subspace $\mathcal{M}$? For example, why is $Y_{t+1}$ found by $\mathcal{P}_{\mathcal{M}}Y_{t+1}$?


title: Why is a projection matrix of an orthogonal projection symmetric?

text: I am quite new to this, so I hope you forgive me if the question is naïve. (Context: I am learning econometrics from Davidson & MacKinnon's book "Econometric Theory and Methods", and they do not seem to explain this; I've also looked at Luenberger's optimization book that deals with projections at an a bit more advanced level, but with no luck). Suppose that I have an orthogonal projection $\mathbb P$ with is associated projection matrix $\bf P$. I am interested in projecting each vector in $\mathbb{R}^n$ into some subspace $A \subset \mathbb{R}^n$. _Question_ : why does it follow that $\bf{P}=P$$^T$, that is, $\bf P$ is symmetric? What textbook could I look at for this result?


title: How to relate the variance of a joint bivariate normal distribution to the variance of a single normal distribution?

text: I have $(X,Y_1)$ with joint bivariate normal distribution. Also, $Y_1$ is conditional on $X$. Therefore $\rho$ is non-zero. Suppose I have $Y_2$ that is also normally distributed and $\mu_{Y_2} = \mu_{Y_1}$ and $\sigma_{Y_2}=\sigma_{Y_1}$. What is the relationship between $\text{Var}(Y_1|X)$ and $\text{Var}(Y_2)$ ? I've been generating samples from $Y_1|X$ and $Y_2$ and it seems that for large sample sizes, the difference $|\text{Var}(Y_1|X)$ - $\text{Var}(Y_2)$| converges to a particular (small) number. This disagrees with theoretical result that I got by calculating $\sigma_{Y_1|X} = \sigma_{Y_1}\sqrt{1-\rho^2}$. What is the interpretation of that?


title: Identifying culs-de-sac (circular housing arrangements) with GIS data

text: I'm thinking about a new project, so I don't have data yet, but I plan on collecting GIS information for houses within a state. Usually in the U.S., these dead-end streets will have a large circle for cars to turn around in, meaning the fronts of each house on the cul-de-sac point toward the center of the circle. In general, my question is how to locate a subset of data clustered in a circular pattern on a 2D grid. If the same technique could allow me to locate dead-ends, that would be useful, but not required. I'll probably need to use a logistic classifier, but I'm not certain what the hypothesis should look like. If this is a well-known problem, I didn't know the keywords to do a proper Google search, so pointers in that direction would be helpful as well.


title: How do I account for within-subject variance to produce the standard error of the mean for different between-subject groups?

text: I have a mixed design with 25 subjects each with 4 repeated measurements in each of three different between subject treatment groups (total n=75). I want to produce the standard error of the mean for each group, but I don't know how to do this in a way that accounts for the within-subject variance. If I simply consider each measurement independent, it artificially decreases the SEM by including too many subjects (i.e. 100 instead of 25), but if I just use 25 it over-estimates the SEM because the variance is considered from 100 "independent" subjects, not 25. I considered simply using the average value for each subject but this feels like "cheating" by simply dismissing this variance. Initially I was using the SE values produced from my SPSS output but these represent an overall SE for the ANOVA values calculated from the model and don't provide the SEM for each group. My colleagues want the group SEM values listed.


title: Is binary hypothesis testing a better statistical term than A/B testing?

text: Is binary hypothesis testing a better statistical term for what business intelligence often refers to as A/B testing? Wikipedia suggests that this is the term used withing the statistics community but Wikipedia is unreliable and there is no tag for such a term here. Since the statistical comparison of multi-grouped randomized experiments predates the term A/B test and the field of business intelligence I would suspect there to be a term. Is it simply "binary hypothesis testing" as suggested by wikipedia? Some historical context would also be appreciated.


title: What do changes in regression coefficients indicate about correlations among predictors?

text: How do you tell if there is a strong, weak, or no correlation between two predictors if you are only given the regression coefficients from two models. One model contains one predictor and the other model contains both predictors.


title: How to handle both text and numbers for PCA in R?

text: I'm relatively new to R and am working with a very large dataset that has a mix of numerical scores (for instance, household income) as well as text values (i.e. race). I was planning on using PCA to analyze this dataset but it only works on numerical data. Is there any good alternative?


title: Modern and practical methods for forecasting time series cross section data

text: I have a data array $X$ which is of size $T \times N \times K$, where $T = 1500$, $N = 1500$ and $K = 10$. Physically, the 1st index $1, 2, \ldots, T$ denotes days, while the 2nd index $1, 2, \ldots, N$ represents locations, and the 3rd index $1, 2, \ldots,K$ represents the $K$ features/variables measured at each location on each day. The dependent variable is another array $Y$ which is of size $T \times N$. Of course, there are correlations in the time and spatial directions. There are also correlations among the $K$ variables/features within each observation. [Some Side Info: For any day $t$, the cross section data can have a multi- level structure, conditional on different grouping criteria. Various grouping criteria can be defined. ] We are interested in forecasting: given all observations until up to $T$, we would like to forecast the Y's on the day $T+1$. Here, some of the $X$ data never vary across time; some other $X$ data do vary across time. Therefore, on day $T$, when we want to forecast $Y$'s on $T+1$, some of the $X$ data $X_{T+1}$ are unknown because they are dynamic and vary across time. What are the modern and practical methods of handling such data? Could anybody please shed some lights on this?


title: How would you explain covariance to someone who understands only the mean?

text: ...assuming that I'm able to augment their knowledge about variance in an intuitive fashion ( Understanding "variance" intuitively ) or by saying: It's the average distance of the data values from the 'mean' - and since variance is in square units, we take the square root to keep the units same and that is called standard deviation. Let's assume this much is articulated and (hopefully) understood by the 'receiver'. Now what is covariance and how would one explain it in simple English without the use of any mathematical terms/formulae? (I.e., intuitive explanation. ;) Please note: I do know the formulae and the math behind the concept. I want to be able to 'explain' the same in an easy to understand fashion, without including the math; i.e., what does 'covariance' even mean?


title: Performing a GLS cross-sectional regression using R

text: Is there a function in R that could perform GLS cross-sectional regression for multiple cases all at once? For example, when regressing stock returns over beta for 100 stocks over a 10-year time series for each case. However, the equation should solve all of the 100 regressions and give only one result, a "multiple regression".


title: Comparing multiple models using the same explanatory variables

text: Using `mgcv` I have made 5 models, each one using a subset of my data, defined by quantiles. The response/explanitory variables are the same in each case, but the deviance explained by each in each model differs, along with the P values. Is there a clever way of comparing the contributions of each explanitory variable in each model? The best I've come up with is using portions of the output of `summary(GAMx)` to create a monstrous table with the `Δ Deviance explained`, `Δ AIC` and `P value` for each each model, with the rows being the smoothed and parametric variables.![monstrous table](http://i.stack.imgur.com/EMG15.png) This is pretty ugly and hard to interpret. Are there any functions or packages in `R` that would allow a comparison of models - with enough detail to show the contributions of each variable? Reproducible code:                library(mgcv)      set.seed(0)      n<-200;sig2<-4      x0 <- runif(n, 0, 1);x1 <- runif(n, 0, 1)      x2 <- runif(n, 0, 1)      y<-x0^2+x1*x2 +runif(n,-0.3,0.3)      g1<-gam(y~s(x0,x1,x2))      g2<-gam(y~s(x0,x1,x2))      g3<-gam(y~s(x0,x1,x2))      g4<-gam(y~s(x0,x1,x2))      g5<-gam(y~s(x0,x1,x2))      All 5 models here will be the same, but that shouldn't actually matter.


title: What code would you put before/after your R session?

text: R allows us to put code to run in the beginning/end of a session. What codes would you suggest putting there? I know of three interesting examples (although I don't have "how to do them" under my fingers here):   1. Saving the session history when closing R.   2. Running a fortune() at the beginning of an R session.   3. I was thinking of having an automated saving of the workspace. But I didn't set on solving the issue of managing space (so there would always be X amount of space used for that backup) Any more ideas? (or how you implement the above ideas) p.s: I am not sure if to put this here or on stackoverflow. But I feel the people here are the right ones to ask.


title: How to perform Prais-Winsten autoregression in SPSS 16?

text: When performing a linear regression on my dataset, Durbin-Watson was very low (0.276). I found a tutorial online that suggested performing an Prais-Winston autocorrelation. The tutoral came with screenshots of SPSS on how to perform the analysis. The screenshots however, are from SPSS 14. I have SPSS 16, and I am not able to find the same test anywhere in my menu structure. Please see the tutorials screenshots (above) and mine (below) here: ![http://i.imgur.com/gN0jG.png](http://i.stack.imgur.com/4ygJ2.png) Can someone direct me to where the Prais-Winsten autoregression can be found in SPSS 16?


title: How to identify models as linear or non-linear?

text: Identify the following models as linear or non-linear. In case of a non-linear model, reduce the model into a linear model by a suitable transformation. $$\eqalign{ (a)\quad&y=\beta_0+\beta_1 x+\beta_2 x^2+e \\\ (b)\quad&y=\frac{x}{\beta_0+\beta_1 x}+e\\\ (c)\quad&y=\frac{\exp(\beta_0+\beta_1 x)}{1+\exp(\beta_0+\beta_1 x)}+e }$$ where $e\sim\mathcal{N}(0,\sigma^2).$ I know that the model of $(a)$ is linear. I think $(b)$ & $(c)$ are non- linear.How can I make them linear models?


title: Get the R2 from GMM Newey-West in SAS

text: I'm running a simple cross-sectional regression where I first run regressions for every year of observations and then I'm running this code to get the Newey-West corrected standard errors:               ods listing close;     ods output parameterestimates=pe;     proc reg data=dset;      by year;      model depvar = indvars; run;     quit;     ods listing;          proc means data=pe mean std t probt;      var estimate; class variable;     run;          proc sort data=pe; by variable; run;          %let lags=3;     ods output parameterestimates=nw;     ods listing close;     proc model data=pe;      by variable;      instruments / intonly;      estimate=a;      fit estimate / gmm kernel=(bart,%eval(&lags+1),0); run;     quit;     ods listing;          proc print data=nw; id variable;      var estimate--df; format estimate stderr 7.4;     run;      What's the best way to get the R2 once I correct for the standard errors for the cross-sectional regression?


title: Can you use composite z-scores in a repeated measures t-test?

text: Can you use composite z-scores to perform a t-test. If so, can you use a composite z-scores in a repeated measures t-test?


title: Marketing Mix Modeling for sales data

text: I have data for 4 years with variables year, month and sales volume only. I want to find the base volume and incremental volume. I want to use Marketing Mix Modeling. Can anybody tell how to define model equation with these variables and how can I use SPSS/R for the same ? Is it possible to build model with only these variables? Is there any other way or statistical method to find these volume ? Please reply... Thanks in advance....


title: What's a good way to use R to make a scatterplot that separates the data by treatment?

text: I'm very new with R and stats in general, but I need to make a scatterplot that I think might be beyond its native capacities. I have a couple of vectors of observations and I want to make a scatterplot with them, and each pair falls into one out of three categories. I would like to make a scatterplot that separates each category, either by colour or by symbol. I think this would be better than generating three different scatterplots. I have another problem with the fact that in each of the categories, there are large clusters at one point, but the clusters are larger in one group than in the other two. Does anyone know a good way to do this? Packages I should install and learn how to use? Anyone done something similar? Thanks


title: Estimating fixed effects for ordinal regression

text: I am running a model when the response variable (behaviour) is an ordinal factor (levels: $0,1,2,3,4$) and I am testing the relationship with two fixed effect factors (year: $A,B,C,D$; colony; $X_1$,$X_2$). I want to extract the estimates of the effect of year and colony on behaviour. In a normal GLM, I understand that this is simple, but for an ordinal model, I am not sure how to incorporate the cut points for the ordinal levels with the fixed effect estimates. From reading about ordinal models, it seems that maybe you get a matrix of probabilities such that you get for colony a $2 \times 5$ matrix and for year a $4 \times 5$. Am I understanding this correctly? And if so, it does not seem possible to extract any kind of residual value for behaviour controlling for fixed effects?


title: Number of needed samples for entropy estimation

text: I would like to estimate the entropy of a source generating binary vectors of length M which are very sparse (just a few 1s), using the naive (empirical) estimator $\hat{H}=-\sum\hat{p}(x)\log\hat{p}(x)$. Is there some result on the number of samples needed to do this with certain precision, or equivalently, the error in calculating this value given finite number of samples?


title: A parellel between LSA and pLSA

text: In the original paper of pLSA the author, Thomas Hoffman, draw a parallel between pLSA and LSA data structures that I would like to discuss with you. **Background:** Taking inspiration the Information Retrieval suppose we have a collection of $N$ documents $$D = \lbrace d_1, d_2, ...., d_N \rbrace$$ and a vocabulary of $M$ terms $$\Omega = \lbrace \omega_1, \omega_2, ..., \omega_M \rbrace$$ A _corpus_ $X$ can be represented by a $N \times M$ matrix of cooccurences. In the _Latent Semantic Analisys_ by _SVD_ the matrix $X$ is factorized in three matrices: $$X = U \Sigma V^T$$ where $\Sigma = diag \lbrace \sigma_1, ..., \sigma_s \rbrace$ and the $\sigma_i$ are the singular values of $X$ and $s$ is the rank of $X$. The LSA approximation of $X$ $$\hat{X} = \hat{U}\hat{\Sigma}\hat{V^T}$$is then computed truncating the three matrices to some level $k < s$, as shown in the picture: ![enter image description here](http://i.stack.imgur.com/snZOE.jpg) In pLSA, choosen a fixed set of topics (latent variables) $Z = \lbrace z_1, z_2, ..., z_Z \rbrace$ the approximation of $X$ is computed as: $$X = [P(d_i | z_k)] \times [diag(P(z_k)] \times [P(f_j|z_k)]^T $$ where the three matrices are the ones that maximize the likelihood of the model. **Actual question:** The author states that these relations subsist:   * $U = [P(d_i | z_k)]$   * $\hat{\Sigma} = [diag(P(z_k)]$   * $V = [P(f_j|z_k)]$ and that the crucial difference between LSA and pLSA is the objective function utilized to determine the optimal decomposition/approximation. I'm not sure he is right, since I think that the two matrices $\hat{X}$ represemt different concepts: in LSA it is an approximation of the number of time a term appears in a document, and in pLSA is the (estimated) probability that a term appear in the document. Can you help me clarify this point? Furthermore, suppose we have computed the two models on a corpus, given a new document $d^*$, in LSA I use to compute it approximation as: $$\hat{d^*} = d^*\times V \times V^T$$   1. Is this always valid?   2. Why I do not get meaningful result applying the same procedure to pLSA? $$\hat{d^*} = d^*\times [P(f_j|z_k)] \times [P(f_j|z_k)]^T$$ Thank you.


title: Example of Stochastic Order

text: I know that stochastic order means the concept of one variable is "bigger" than another. But I would like some example to visualize/understand it better.


title: What are the optimizations or goals to consider when using stratified sampling?

text: In stratified sampling, what are the optimization considerations? For example, the sample size per stratum could be defined with proportional allocation as $n_h=n\frac{N_h}{N}$, where $N_h$ is the population size for stratum $h$, $N$ is the total population size, $n$ is the total sample size, and $n_h$ is the sample size for stratum $h$. For a binary variable, I can choose $n$ so that my margin of error $e$ is 1%, 5% or 10%. $$e=z_{1-\frac{\alpha}{2}} \sqrt{\frac{pq}{n}}$$ I could also choose my $n = \sum_{h=1}^H n_h$ such that the margin of error for each stratum $e_h$ is 1%, 5% or 10%. $$e_h=z_{1-\frac{\alpha}{2}} \sqrt{\frac{pq}{n_h}}$$ But my question/concern is, how many different considerations are there and which ones are the most important ones in the context of stratified sampling?


title: Why is r-squared so small when there does seem to be a dependency in my data?

text: What I need to do is to find a model that can predict what the **observation** should look like given a **factor** input. I am doing a simple linear fit in R (i.e., `lm(observation~0+factor, data=d)`); the $R^2= 0.002$, which is really small. However, when I do a `SELECT AVG observation by 0.001 BRACKET factor`, the result is something like:               factor   | average observation     ----------------------------------------     -0.003        -2     -0.002        -2     -0.001        -1      0.000         1      0.001         0      0.002         1      0.003         2      It definitely seems to me that there is a pattern here, but somehow this pattern is not captured by a linear model. Is my understanding correct?


title: How to report a SVM model to a 3rd party after cross-validation?

text: I have a binary classification problem. I trained my dataset using a Support Vector Machine (SVM). Now I want to report the model I trained to a 3rd party so that they can use. For the primal probem of SVM, the classifier equation is of the form,               y = w.x + b      where `w` is the weight vector, `b` is the constant and `x` is the new instance. Depending on the value of `y`, we can predict the class. I could report the weight vector w and constant b which was obtained my model. But if I want to test my model using cross-validation, then how do I report my model? Everytime I cross-validate, the equation of the classifier built will change.


title: Parameter Estimation of a Poisson mixture model

text: I want to estimate the parameters for a Poisson mixture model with 2 (and later 3) Poisson distributions. I want to use Matlab and have numerical problems to solve the loglikelihood of the mixture model (values get to high or low). Does anyone know how to handle this problem? I already tried to write a similar function like logLikePoisMix (https://r-forge.r-project.org/scm/viewvc.php/pkg/HTSCluster/R/logLikePoisMix.R?view=markup&root=htsfilter&sortby=rev&pathrev=60), which uses a "calculation trick from Panos" (?), but this also didn't work. Thanks for your help!


title: Optimizing False Negative Rate after Logistic Regression

text: I created a probit model and tested it against a random sub sample of my dataset. I am interested specifically in seeing how many data points I can predict to be FALSE without having too many that are actually TRUE. Using the threshold of 0.1 (see below), I was able to predict to be FALSE about 30% with a false negative rate of 2.5%.   However, I don't know if this is optimal. Is there a way for me to pick my threshold that maximizes my FALSE predictions while minimizing my false negatives?               glm.fit=glm(Outcome~A+B+C+D+E+F+G,data=myData,family=binomial(link="probit"))     test=mysample <- myData[sample(1:nrow(myData),10000,replace=FALSE),]     glm.probs =predict(glm.fit,test, type="response")     glm.pred=rep(0,10000)     glm.pred[glm.probs>.1]=1     x <- sum(glm.pred == 0 & test$Outcome == 1)     y <- sum(glm.pred == 0)


title: Probability density function within [0,1] with specifiable mode

text: I needed a probability density function which worked on the interval $[0,1]$, had kind of a bell shape, and had an adjustable mode / peak $p$. I thought of a pdf $f(x|p)$, given by \begin{equation} f(x|p) = \left\\{ \begin{array}{l l l } \frac{ (x^{- \ln 2/\ln p})^2 \cdot (1-x^{- \ln 2/\ln p})^2 }{ \log(p) \left( \frac{x^{1-\frac{4 \log 2}{\log p}}}{\log(p) - 4\log(2)} -\frac{2 x^{1-\frac{3 \log 2}{\log p}}}{\log(p) - 3\log(2)} +\frac{x^{1-\frac{2 \log 2}{\log p}}}{\log(p) - 2\log(2)} \right) } &\quad \text{ for } 0<x<1 \\\ 0 &\quad \text{ otherwise} \end{array} \right. \end{equation} which   * has a peak at $x=p$ for $0<p<1$.   * $P(X\le 0) = P(X\ge 1) = 0$   * has a shape similar to the bell shape   * looks skewed to the left for $p>1$ and to the right for $0<p<1$ Or equivalently: \begin{equation} f(x|p) = \left\\{ \begin{array}{l l l } \frac{ (x^a)^2 \cdot (1-x^a)^2 } { (4a+1)^{-1} - 2(3a+1)^{-1} + (2a+1)^{-1} } &\quad \text{ for } 0<x<1 \\\ 0 &\quad \text{ otherwise} \end{array} \right. \end{equation} which which has its peak at $p=-\frac{\log 2}{\log x}$ Is there a similar pdf (or exactly this one) used in literature? What is it called? **PS** : note that the given pdf is not symmetric: $f(x|1-p) \neq f(1-x |p) $


title: R as a backend processing engine

text: I have an R script that process a CSV file and spit out result as CSV. This work fine in desktop environment. Now I would like to implement some think like this:   1. set up an email where user could send CSV file as attachment   2. On a server, have an program to check email automatically and pass the CSV file to a specific folder (I don't know exactly how to do this, but I have people in the team to do this)   3. Every time when there are new CSV file come into the folder, get R to pick up this CSV file and run the script. The result is a new CSV file.   4. Deliver the result back to the user by sending an email. I'm kind of stuck on number 3. Did anyone setup R (or R server) to run in this fashion. Any suggestion is much appreciated.


title: predict function and categorical variables in R

text: This is more of a general question about how the predict function treats categorical variables and how to interpret the output from predict. I have a zeroinfl model to predict the number of animals encountered:               b9 <- zeroinfl(Count ~ as.factor(Area) + as.factor(Season) | 1, dist="negbin", data=total)      where Count is the number of animals and the explanatory variables are Area and Season, both are coded as factors in the model. Area has three levels and Season has 4 levels. Coding the two as factors allows for R to create dummy variables for each variable for use in the model. When I use the predict function to predict the number of animals for a larger data set, I want to make sure I'm understanding what is happening. My newdata for predict is: `newdata<- as.data.frame(Season, Area)`. Both variables are coded as factors and the dataframe is in a long format. There are records for each combination of Season and Area that correspond to trips taken over the course of 8 years. There are 113,804 rows of data in the newdata data,frame.               str(newdata)     'data.frame':   113804 obs. of  2 variables:     $ Season: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...     $ Area  : Factor w/ 3 levels "625","631","Bay": 3 3 3 3 3 3 3 3 3 3 ...      Example:                  Season Area     1       1  Bay     2       1  625     3       1  631     4       2  Bay     5       2  625     6       2  631     7       3  Bay     8       3  625     9       3  631     10      4  Bay     11      4  625     12      4  631        1. Do I need to create dummy variables for all levels for the two variables for input into predict, or does predict function act like zeroinfl where if the variables are coded as factors, this is done automatically by R.   2. Predict returns values of 0.0461 - 0.6015. If I am trying to predict the number of animals how do I interpret this? Since no predicted values are greater than 1 and I need whole numbers, I rounded the predicted data so that any value less than 0.5 was equal to 0 and any value greater than 0.5 was equal to 1. Does this seem correct? Thanks for any help.


title: Meaning of "design" in design matrix?

text: In linear regression, $Y= X\beta$, why is $X$ called the design matrix? Can $X$ be designed or constructed arbitrarily to some degree as in art?


title: Computational Complexity of Prediction using SVM and NN?

text: I've seen answers discussing the complexity of training SVMs and neural nets, but how about for _predicting new responses_ once a model has been trained? For context, I'm working on an app that should produce predictions in near real-time given incoming pixel data and I'm looking for a machine learning algorithm that can handle complex separating planes and predict as fast as possible.


title: What are the units in this biplot?

text: This is a plot of my data ![enter image description here](http://i.stack.imgur.com/OaPQx.png) These are the values:                  xvalues  yvalues     1   1.091186     2   2.653722     3   3.309146     4   5.206479     5   5.115582     6   8.537005     7   10.013147     8   9.802291     9   10.667769     10  5.809750     11  9.624475     12  11.806013     13  13.587066     14  14.146781     15  13.707472     16  12.891355     17  19.435301     18  16.122108     19  17.768536     20  23.813027     21  21.819081     22  23.556074     23  21.170983     24  27.621148     25  22.932580     26  20.704689     27  25.530339     28  26.227371     29  26.051016     30  31.047145      I do now a PCA and a biplot of it ![enter image description here](http://i.stack.imgur.com/IqWXm.png) According to Jeromy Anglim in: Interpretation of biplots in principal components analysis in R _The left and bottom axes are showing the loadings; the top and right axes are showing principal component scores._ I want to be sure I really get this. Lets start with the loadings: these can be visualized in R by writing:               results <- prcomp(your_data)     results$rotation                     PC1        PC2     xvalues   0.7235616  -0.6902599     yvalues   0.6902599   0.7235616          summary(results)          Importance of components:                            PC1     PC2     Standard deviation     12.0747 1.56606     Proportion of Variance  0.9835 0.01654     Cumulative Proportion   0.9835 1.00000      Now lets look at the red arrow of xvalues. Its tip is around 0.25 in the x-axis of the loadings. But according to the loadings I have just writen, it should be around 0.72. What am I missing? Finally, lets look at the point 1. According to the axes, that is telling me the principal component score. Is it that the coordinate in the new frame of reference? It doesn't make sense to me because I think that the new origin of the axes is around the point (15,15) in the Plot 3. If I look at that (and I guess that I am completely wrong here), the point one should have a coordinate around -20 or so, and not 40. Where is my mistake? **Update** I tried plotting this:               plot(pca_results$x)      ![enter image description here](http://i.stack.imgur.com/eyNwX.png) Here it can be seen that the first point has the coordinate that I thought it had to have. But, still, what are the units in the biplot then??


title: Stata rreg and reg, robust

text: I have a general question about regression in Stata. What is the main difference between commands: `rreg` (Huber regression) and `reg, robust`? What is the explanation for one and the other?


title: What is filtration?

text: I am studying Brownian motion and came across the concept of Filtration. However I can't understand how this concept relates to Brownian motion. My notes contained some gibberish about being measurable and adpated. But all these concepts are new to me that I'm having a hard time to grasp what they are about. Can somebody clear that up for me?You should assume that I know nothing about these concepts as I'm totally lost.


title: Any suggestions for a good undergraduate introductory textbook to statistics?

text: I am hoping you can give me some suggestions. I am teaching in a very diverse (made of minority groups) college and the students are mostly Psychology majors. Most students are fresh from high school but some of them are older returning students above 40. Most of the students have motivational problems and aversion to math. But I am still looking for a book that covers the basic curriculum: from descriptive to sampling and testing all the way to ANOVA, and all in the context of experimental methods. The department requires me to use SPSS in class, but I like the idea of building the analysis in a spreadsheet such as excel. p.s. the other teachers use a book that I don't like because of the extensive reliance on computational formulae. I find using these computational formulas - rather than the more intuitive and computationally intensive formula that is consistent with the rational and basic algorithm- unintuitive, unnecessary and confusing. This is the book I refer to Essentials of Statistics for the Behavioral Sciences, 7th Edition Frederick J Gravetter State University of New York, Brockport Larry B. Wallnau State University of New York, Brockport ISBN-10: 049581220X Thank you for reading!


title: Comparing two means

text: I have a question on comparing means and I would be really thankful if you could help... I have a data set that contains responses from a questionnaire (so one group only, N=216). The questionnaire forms 3 scales (A,B,C). My means for these scales are: A:mean=14.92, with SD=6.96 B:mean=10.51, with SD=7.32 C:mean=14.04, with SD=11.12 I would like to compare these three means to see which scale is higher (where people scored higher). What should I do? Should I first transform my data into z-scores before making comparison? Thank you for your time and help.


title: ANCOVA in observational studies: what are the assumptions?

text: Using ANCOVA when groups differ on the covariate is controversial, although Tabachnick and Fidell write that this is a plausible function of ANCOVA in quasi-experimental (or observational) studies. As they state: > The second use of ANCOVA commonly occurs in nonexperimental situations when > subjects cannot be randomly assigned to treatments. ANCOVA is used as a > statistical matching procedure, although interpretation is fraught with > difficulty [...]. ANCOVA is used primarily to adjust goup means to what they > would be if all subjects scored identically on the CV(s). Differences > between subjects on CVs are removed so that, presumably, the only > differences that remain are related to the effects of the grouping IV(s). > (Differences could also, of course, be due to attributes that have not been > used as CVs.) This second application of ANCOVA is primarily for descriptive > model building: the CV enhances prediction of the DV, but there is no > implication of causality. If the research question to be answered involves > causality, ANCOVA is no substitute for running an experiment. Moreover, in this question the same issue was addressed, and the use of ANCOVA for intact groups was encouraged. **My question is: in these situations, in which the assumption of independence of the covariate from the treatment variable is violated, what are the assumptions? For example, must the covariate be correlated with the dependent variable inside the groups? Or are the assumptions simply the same as for ANOVA?**


title: How to estimate the probability of a rare event about which observations can only be made in quantized time?

text: I have reduced the estimation problem of a real event (a technical failure happening, the fact of which is checked in regular time intervals) to the following problem: we have a non-fair coin, which gives a head for almost every throw. If gives a tail extremely rarely. Let's say we throw the coin every second. After an hour, a tail comes up. After two minutes, another. After 2 hours, another. How can we estimate the probability of a tail occurring, and how can we have a good guess at the reliability of the estimate after a given number of (e.g. few dozen) occurrences of a tail? My problem is that the event is rare enough that it's very hard to get a reliable measurement for some of the longest time periods when it doesn't occur (it just takes a lot of time).


title: The best distribution and regression method when modelling count data?

text: I am trying to fit a regression model where the dependent variable is number of messages received ( **MsgReceived** in sample data below) by an individual/user, and the independent variables are a mix of count and binary variables like " **TimeActive** ", " **BioAvailable** " etc. Here, **TimeActive** is a count of weeks the user was active and **BioAvailable** is a 1-0 identifier stating whether the user has filled out his Bio Here's some example data:               TimeActive,BioAvailable,X1,X2,X3,X4,X5,X7,X8,X9,X10,X11,Count1,Count2,Count3,MsgReceived          35,1,1,0,0,0,1,1,0,1,0,0,3,0,3,16       34,1,1,0,1,1,0,1,0,1,0,0,20,23,37,11       34,0,0,0,1,0,0,1,0,1,1,1,6,8,22,19       35,0,0,0,1,1,1,1,0,1,0,0,3,23,5,13       32,0,0,0,1,0,1,1,0,1,1,1,0,75,11,40       0,0,0,0,0,0,0,1,0,1,1,1,0,0,0,7       21,0,0,0,0,0,0,0,0,0,0,0,3,33,39,97       13,1,1,0,0,0,0,1,0,0,1,1,1,0,0,12       34,0,0,0,1,1,0,1,0,1,0,0,35,52,2,37       33,1,0,0,1,1,0,1,0,1,0,0,0,9,16,136       31,1,1,0,1,1,0,1,1,1,1,0,5,1,12,46       0,0,0,0,1,0,1,1,0,1,1,1,0,3,8,20       29,1,1,0,1,1,1,1,0,1,0,0,44,161,45,8        I am wondering if fitting a generalized linear model using a Poisson distribution is still the best fit even though the count of messages is over the life of a user's activity and not just a session. Assuming that a user's lifetime is a period seems fair. Is this correct? If not, what distribution and regression method would be a better fit?


title: Interaction term significant but in the plot the lines don't intersect

text: I tested the moderating effect of the variable X' on the relation between X and Y. The results were positive and significant. However, when I plotted the results, I found out the two lines become very, very close at the higher levels, but do not "cut each other". I wonder if these results make sense... I'm trying to clarify if X' in fact moderates the relation between X-Y (the plot is confusing me). Thanks!


title: Covariance matrix in Gaussian Discriminant Analysis

text: I am at a point where I stumbled upon something I don't fully understand. I have to compute the covariance matrix using this formula: $$\Sigma=\frac{1}{m}\sum_{i=1}^{m}(x^{(i)}-\mu_{y^{(i)}})(x^{(i)}-\mu_{y^{(i)}})^{T}$$ with $$\mu_0=\frac{\sum_{i=1}^{m}1\\{y^{(i)}=0\\}x^{(i)}}{\sum_{i=1}^{m}1\\{y^{(i)}=0\\}}$$ and $$\mu_1=\frac{\sum_{i=1}^{m}1\\{y^{(i)}=1\\}x^{(i)}}{\sum_{i=1}^{m}1\\{y^{(i)}=1\\}}.$$ My question, $\Sigma$ is a matrix right? $\mu_0$ and $\mu_1$ are also matrices. And here comes what I don't get. $x^{(i)}$ is a training example, let's say with 2 features, say width and height. $\Sigma$ should be a square matrix.


title: How to quantify the significance of the difference between two z-scores?

text: I have one sample and several features. I calculate a z-score for various features, and for various combinations of features. Is there a way to quantify the significance of the difference between two z-scores? For example, suppose I calculate a z-score for the combination of features (male, homeowner) and also for the combination (male, homeowner, married). For one scenario, suppose the first z-score is 3.0, and the second z-score is 6.0. For a different scenario, suppose the first z-score is 25.000, and the second z-score is 25.001. I'm trying to determine whether or not (male, homeowner, married) is a better predictor than (male, homeowner). Is there a generally applicable formula to quantify the significance of the difference, perhaps taking into account the sample size? I have no training in statistics, and am trying to implement a software program. A cookbook formula accompanied by some beginner-appropriate explanation would be ideal and greatly appreciated. Thanks


title: Using the R forecast package with missing values and/or irregular time series

text: I am impressed by the R `forecast` package, as well as e.g. the `zoo` package for irregular time series and interpolation of missing values. My application is in the area of call center traffic forecasting, so data on weekends is (nearly) always missing, which can be nicely handled by `zoo`. Also, some discrete points may be missing, I just use R's `NA` for that. The thing is: all the nice magic of the forecast package, such as `eta()`, `auto.arima()` etc, seem to expect plain `ts` objects, i.e. equispaced time series not containing any missing data. I think real world applications for equispaced-only time series are definitely existent, but - to my opinion - v e r y limited. The problem of a few discrete `NA` values can easily be solved by using any of the offered interpolation functions in `zoo` as well as by `forecast::interp`. After that, I run the forecast. My questions:   1. Does anyone suggest a better solution?   2. **(my main question)** At least in my application domain, call center traffic forecasting (and as far as I can imagine most other problem domains), time series are not equispaced. At least we have recurring "business days" scheme or something. What's the best way to handle that and still use all the cool magic of the forecast package?  Should I just "compress" the time series to fill the weekends, do the forecast, and then "inflate" the data again to re-insert NA values in the weekends? (That would be a shame, I think?) Are there any plans to make the forecast package fully compatible with irregular time series packages like zoo or its? If yes, when and if no, why not? I'm quite new to forecasting (and statistics in general), so I might overlook something important.


title: Comparing regression coefficients of models with differnet dependent variables

text: I am interested in comparing response (dependent) variables (2+) in a regression format where the models share the same data and design matrix. As an example, lets say I am interested in whether numerical counts of aquatic organisms, or the biomass of those organisms is more "sensitive" to the same disturbance gradient (e.x., flood frequency). So I figured, comparing regression coefficients, specifically slopes would be appropriate. Could you please tell me if you have a better idea, or if this is not appropriate.   Thank you.


title: Is there a term to indicate data outside the control/specification limits?

text: Is there a term that indicates a data that has fallen outside CL/SL for Statistical Process Control? Most seem to call it simply "Data out of bounds of CL/SL" but I'd like to know if there's a specific term for it.


title: Pearson residuals formula in a multinomial logit model

text: Could someone tell me what's the formula for the Pearson residuals in a multinomial logit model? I tried to look for it but I haven't found anything.


title: How do I show that $(\hat{\beta_1}-\beta_1)$ and $ \bar{u}$ are uncorrelated, i.e. that $E[(\hat{\beta_1}-\beta_1) \bar{u}] = 0$?

text: I asked this question and thought I found a solution to it by working through the steps in my textbook. Unfortunately the comments told me that I was completely wrong and only got the proof to work by chance, so now I'm really lost and frustrated. To fix notation, the model is $$ y_i = \beta_0 + \beta_1 x_i + u_i, $$ where $u_i$ is the error term. Earlier I showed that $\hat{\beta}_1 = \beta_1 + \sum_{i=1}^n w_i u_i$, where $w_i = \frac{x_i - \bar{x}}{SST_x}$, and $SST_x = \sum_{i=1}^n (x_i - \bar{x})^2$. Here is what I did so far: \begin{align} E[(\hat{\beta_1}-\beta_1) \bar{u}] &= E[\bar{u}\displaystyle\sum\limits_{i=1}^n w_i u_i] \\\ &=\displaystyle\sum\limits_{i=1}^n E[w_i \bar{u} u_i] \\\ &=\displaystyle\sum\limits_{i=1}^n w_i E[\bar{u} u_i] \\\ &= \displaystyle\sum\limits_{i=1}^n w_i \left(Cov(\bar{u}, u_i) + E(\bar{u})E(u_i)\right) \\\ &= \displaystyle\sum\limits_{i=1}^n w_i Cov(\bar{u}, u_i) \\\ &= \displaystyle\sum\limits_{i=1}^n w_i Cov\left(\frac{1}{n} \displaystyle\sum\limits_{i=1}^n u_i, u_i\right) \\\ &= \frac{1}{n}\displaystyle\sum\limits_{i=1}^n w_i Cov\left(\displaystyle\sum\limits_{i=1}^n u_i, u_i\right) \\\ &= \frac{1}{n}\displaystyle\sum\limits_{i=1}^n w_i E\left(\left[u_i - E(u_i)\right]\left[ \displaystyle\sum\limits_{i=1}^n u_i - E\left(\displaystyle\sum\limits_{i=1}^n u_i\right)\right]\right) \\\ &= \frac{1}{n}\displaystyle\sum\limits_{i=1}^n w_i E\left(\left[u_i - E(u_i)\right]\left[ \displaystyle\sum\limits_{i=1}^n u_i - \displaystyle\sum\limits_{i=1}^n E\left(u_i\right)\right]\right) \\\ &= \frac{1}{n}\displaystyle\sum\limits_{i=1}^n w_i E\left(\left[u_i - 0\right]\left[ \displaystyle\sum\limits_{i=1}^n u_i - 0\right]\right) \\\ &= \frac{1}{n}\displaystyle\sum\limits_{i=1}^n w_i E\left(u_i\displaystyle\sum\limits_{i=1}^n u_i\right) \\\ &= \frac{1}{n}\displaystyle\sum\limits_{i=1}^n w_i \left[E\left(u_i u_1 +\cdots + u_i u_n \right)\right] \\\ &= \frac{1}{n}\displaystyle\sum\limits_{i=1}^n w_i \left[E\left(u_i u_1\right) +\cdots + E\left(u_i u_n \right)\right] \\\ \end{align} At this point, I want to say that because we assume the errors for each observation are i.i.d: \begin{align} &= \frac{1}{n}\displaystyle\sum\limits_{i=1}^n w_i \left[E(u_i) E(u_1) +\cdots + E(u_i) E(u_n)\right] \\\ \end{align} But if this is the case then since $E(u_i) = 0, \forall i$, everything just cancels and I don't need to use the definition of $w_i$, which one of the comments mentioned I needed to use last time. This is in chapter 2 of my book and the book says that "it's obvious that..." so maybe it is obvious. I keep asking for hints but I feel like I don't understand them, but I want to keep trying. Is this correct? If not, does anyone have a hint on this part of the problem, assuming I'm on the right path? I can probably make the steps less verbose, but I want to make sure it's right first. The last time I worked this out, I reached the right conclusion but the proof was all wrong, so I'm assuming I only have it right by chance again.


title: Combining two unequal normal distributions

text: Let $X_1$ and $X_2$ be independent, normal distributed random variables with equal mean $\mu$ but non-equal standard deviations $\sigma_1$ and $\sigma_2$. Suppose I know $\sigma_1$ and $\sigma_2$ and I have $n$ samples $x_{11},\ldots,x_{1n}$ from $X_1$ and m samples $x_{21}, \ldots, x_{2m}$ from $X_2$, what is the best estimator for $\mu$? What is its distribution? (edit: I'm mainly interested in the n=1, m=1 case)


title: How to get the p-value for the full model from R's coxph?

text: In a coxph model with 2 variables (let's say age and sex), how can I get the p-values for LR, Wald and score tests which I see in the summary(coxphobject) as below?               Concordance= 0.653  (se = 0.058 )     Rsquare= 0.085   (max possible= 0.924 )     Likelihood ratio test= 30.09  on 12 df,   p=0.002708     Wald test            = 34.73  on 12 df,   p=0.0005169     Score (logrank) test = 34  on 12 df,   p=0.0006736,   Robust = 17.61  p=0.1279      I think those values represented in the summary as `p=...` for each model are the p-values for the full model, but I cannot figure out how to retrieve them. If there is only one variate, let's say age, then I can retrieve the p-value using `summary(coxphobject)$coefficients[5]`, which is the same as the one I see for the Wald test in the `summary(coxphobject)`. But when there are two variables (age and sex), `summary(coxphobject)$coefficients[1,5]` and `summary(coxphobject)$coefficients[2,5]` give the p-values for age and sex separately, but they both are different from the p-values I see for the full model in the `summary(coxphobject)` for LR, Wald and score tests. I appreciate any help. Thanks.


title: What exactly is the vector fitting on a nMDS plot telling me about community structure?

text: I think I have a basic misunderstanding of what the vectors typically plotted over an nMDS (non Metric Multidimensional Scaling) plot are telling me. I understand the direction of the vectors point to the direction of greatest change in the plot. But it was my understanding the length of each arrow (calculated from the R2 values) is based on the importance each vector has on the makeup of the plot. Thus a long arrow has a large effect on the position of each point in the plot and vice versa. But here is my concern- the following made up data represents 8 species at 10 different sites. The 8 species represent different extremes of possible community change.               library(vegan)           # High Abundance, Low Change     HALC <- c(100,99,101,98,99,100,100,101,99,100)      # Low Anundance, Low Change     LALC <- c(1,2,1,2,2,2,1,2,1,2)      # High Abundance, High Absolute Change     HAHAC <- c(100,99,98,91,86,50,45,30,21,9)      # High Abundance, High Absolute Change     LAHAC <- c(9,21,30,45,50,86,91,98,99,100)      # High Abundance, Changes 10%     HA10C <- c(100,101,102,103,104,104,106,107,108,110)      # Low Abundance, Changes 10%     LA10C <- c(1,1.01,1.02,1.03,1.04,1.04,1.06,1.07,1.08,1.1)      # High Abundance, Changes 100%     HA100C <- c(100,110,120,130,140,140,170,175,187,200)      # Low Abundance, Changes 100%     LA100C <- c(1.00,1.10,1.20,1.30,1.40,1.40,1.70,1.75,1.87,2.00)          DF <- c(HALC,LALC,HAHAC,LAHAC,HA10C,LA10C,HA100C,LA100C)     DF <- data.frame(HALC,LALC,HAHAC,LAHAC,HA10C,LA10C,HA100C,LA100C)     row.names(DF) <- c(1,2,3,4,5,6,7,8,9,10)           # Not sure what is a better transformation,'total' or 'norm'     RWideHab <- decostand(DF, method='norm')          set.seed(999)     mdsHab <- metaMDS(RWideHab, dist='bray',trymax=99)           # A function to calculate importance of each species using Procustes analysis       #  (from Jari Oaksanen in a different question)     impo <- numeric(8)     for(i in 1:8) impo[i] <- procrustes(mdsHab, metaMDS(RWideHab[,-i]))$ss     IMPO <- data.frame(impo)     row.names(IMPO) <- c('HALC','LALC','HAHAC','LAHAC','HA10C','LA10C','HA100C','LA100C')           # plot NMDS as normal     plot(mdsHab, display='sites',type='t')      # Use envfit to get environmental vectors     envfitHab <- envfit(mdsHab, RWideHab,perm=999)      # overlay arrows on above plot     plot(envfitHab)          VectHab  <- as.data.frame(scores(envfitHab,display='vectors'))     AbVecHab <- abs(VectHab)     Sweep    <- sweep(AbVecHab, 2, colSums(AbVecHab), "/")     DFVec    <- data.frame(envfitHab$vectors$arrows, envfitHab$vectors$r,                             envfitHab$vectors$pval)     Env_Vector <- data.frame(DFVec,VectHab,Sweep)          colnames(Env_Vector)<- c("NMDS1_Raw", "NMDS2_Raw", "R2", "p", "NMDS1_Scaled",                               "NMDS2_Scaled", "NMDS1_Relative", "NMDS2_Relative"  )     IMPO     Env_Vector      The variable `IMPO` should be a measure of importance calculated from Procrustes analysis, so the higher the Procrustes number, the more important that species is in determining the make up of the plot. Yet in this case, the Procrustes numbers don't tend to agree with the vector lengths in the plot (vector length data is found in the variable `Env_Vector`) . Especially with the species designated as `LALC` (Low Abundance, Low Change), which ranks high-ish in IMPO, but has a very short vector. So the `IMPO` (Procrustes analysis) and vector fitting don't seem to agree, and I am assuming I am misunderstanding what exactly the vector magnitude information is telling me about the importance of different species in my plot.


title: How to reload the WinBUGS computation results to bugs object again?

text: I have run a very long computation in WinBUGS (million iterations) using the R2WinBUGS package from within R:               bugs.object <- bugs(...)      but the R crashed. **How do I reload the bugs.object into R again without running winbugs again?** I tried this (I have 3 chains):               out <- read.bugs(paste("coda", 1:3, ".txt", sep = ""))      but the `out` data structure is completely different from the bugs object (as it is, it is unusable). I tried to convert it with `as.bugs.array`:               bugs.object <- as.bugs.array(out, model.file = "ttest.txt", n.iter = 1000000, n.burnin = 300000, n.thin = 2, program = "WinBUGS")      but it doesn't work.


title: How can this series pass the structural change test?

text: I have this series:               r <- c(14.1649994,17.3997831,32.8877030,6.1547551,6.3594748,101.9584077,196.8411517,23.8292043,11.7583992,2.2915478,5.5088687,16.8188390,2.3077979,2.6666377,2.0955362,0.9319688,0.6688078,1.0523184,6.8158264,13.1813211,17.8467944,1.9876636,7.0219306,3.9272217,4.3325448,7.8778125,1.1878058,1.7207762,2.9655665,5.8571428)      I'm using the `strucchange` package to detect the structural changes. The result is:               > sctest(r~1)              Recursive CUSUM test           data:  r ~ 1       S = 0.7785, p-value = 0.1579      How is it possible that the test didn't detect changes? There is a very big move rising around 200. Take a look at the plot of the series below. ![enter image description here](http://i.stack.imgur.com/iwvB2.png)


title: Scaling the backward variable in HMM Baum-Welch

text: I am just trying to implement the scaled Baum-Welch algorithm and I have run into a problem where my backward variables, after scaling, are over the value of 1. Is this normal? After all, probabilities shouldn't be over 1. I am using the scale factor I obtained from the forward variables: $$ c_t = 1 / \sum_{s\in S}\alpha_t(s)\\\ $$ where c_t is the scaling factor for time t, alpha is the forward variable, s are the states in the hmm. For the backward algorithm I implemented it in java below:               public double[][] backwardAlgo(){             int time = eSequence.size();             double beta[][] = new double[2][time];                  // Intialize beta for current time             for(int i = 0; i < 2; i++){                 beta[i][time-1] = scaler[time-1];             }                  // Use recursive method to calculate beta             double tempBeta = 0;             for(int t = time-2; t >= 0; t--){                 for(int i = 0; i < 2; i++){                     for(int j = 0; j < 2; j++){                         tempBeta = tempBeta + (stateTransitionMatrix[i][j] * emissionMatrix[j][eSequence.get(t+1)] * beta[j][t+1]);                     }                     beta[i][t] = tempBeta;                     beta[i][t] = scaler[t] * beta[i][t];                     tempBeta = 0;                 }             }             return beta;         }      The scales are stored in the array called scaler. There are 2 states in this hmm. I should also note that the scale factors I am getting are over 1 as well.


title: Poisson random variable conditioned on sum

text: If $X_1,X_2...X_N$ are independent Poisson variables with parameters $\lambda_1,\lambda_2...\lambda_3$, then given $\sum_iX_i=N$, we have that $X_i \sim \mathrm{Binom}(N,\frac{\lambda_i}{\sum_j\lambda_j})$ (straight from wikipedia!) I am looking for a proof of this -- I haven't really studied the relevant stats to know how to attack this problem. Any help would be great. Cheers!


title: Linear Mixed Model Construction Suggestion

text: I'm currently working on some data that requires me to use linear mixed model. **Explanation of Data**   * Experiment of Drugs on Mouse   * I have 3 drugs and 1 control   * Unbalanced number of mouse for each drug and control (4~13)   * Measured by Day (range from 0 to max of 33)   * The measured aspect is Tumor Volume **My thoughts**   * Fixed effect: Day*Drug interaction   * Random effect: ~Day|MouseID (MouseID is of course the random intercept since each mouse have different tumor volume; and it is reasonable to think that each tomor has a unique growth rate on each mouse)   * Check the output:      * Random effect StdDev [Can I calculate how much of the residual is explained by the random effect -> (Intercept+Slope) / (Intercept+Slope+Residual)]     * Fixed effect: check p-values for drug and interaction     * Correlation: I don't know what can we know from this part **My Question**   * Is it reasonable to calculate how much of the residual is explained by the random effect?   * How to interpret correlation part?   * Are there any other factors or situations that I need to think of?   * Which number or part is the most important in the summary of a lme model, and how should we interpret that?    * Are there any better model construction ideas? Thanks in advance!!


title: Least square regression with L1 regularization and non-negativity constraint

text: There are two functions associated by the model > $a(x) = \int_{k_1}^{k_2} b(k)\exp(-kx)dk$ where $a(x)$ is the experimental data I have, and $b(k)$ is the information I want to get. Or I can write in a matrix form if k space is discretized: > $a = Mb$ Now I would like to solve this inverse problem with L1 regularization with b>=0 constrain: > min: $|Mb-a|^2 + \lambda|b|$ > > subject to: $b\geq 0$ with $\lambda$ as the regularization parameter. (Actually this problem is similar to a previous one, but this time it is L1 norm regularization.) My questions are:   1. While general solvers for the unconstrained least square with L1 regularization problems exist, I have some difficulties finding a solver for the constrained case like here. Is there an existing solver/routine I can make use of? If not, is there an easy way to get on this problem?    2. I don't know the noise level. Someone suggested cross-validation method to get the regularization parameter. I have some difficulties understanding this method. Could anyone give an easy explanation and point me to the appropriate resources for further understanding?  Thank you so much! Update: regarding the first question above, I have found this Matlab solver is really great.


title: Regression of irregular time series

text: I am beginning to look at high frequency data sets. I have only ever dealt with regularly spaced data up until now so wondering what best practice is. If I have two time series X and Y with Y being the dependent variable, if both time series are irregular is the only way to regress Y against X to perform some form of interpolation on both series? Or is there a way to consider the time between data points as equal in duration and then perform a regression with X lagged against Y. I assume would have to consider that if you get a data point for X at time stamp 0 mics then the next data point in Y must occur after xx mics since it takes time for any response to occur in Y in relation to X. Is there a good way to determine that duration of response time? or is this just an assumption?


title: What does more than chance mean?

text: I'm interested in what exactly sentences like this mean: > In the former case, participants correctly chose the larger city in 55% of > the pairs, in the latter in 54% of the pairs, both rates significantly > different from chance. Is it basically just that its significantly different from 50%? UPDATE: Here's another one:               But now here's another one I've found: In a pretest, one brand had been rated as higher quality, and participants could identify the higher-quality product 59 percent of the time in a blind test (substantially higher than chance, which was 33 percent). I still don't get how is 33 determined?      How do they get 33?


title: Separating two complex-valued datasets that have been multiplied together

text: I have two complex-valued datasets, **A** and **B** that can be considered as vectors with the same number of elements. The datasets are multiplied together using complex point-by-point multiplication, such that the first element of **A** is multiplied by the first element of **B** , the second element of **A** is multiplied by the second element of **B** , and so on. The multiplication is similar to the Matlab dot operator, `C = A.*B`. The resulting dataset is **C**. **A** is a filter kernel in the frequency domain, which models the effects of a medium on a signal **B**. **B** is the reflection of the signal from the medium. I want to treat **B** as noise. The point-by-point multiplication of **A** by **B** in the frequency domain is equivalent to a convolution in the time domain. Suppose that I do not know **A** or **B** , but using a statistical algorithm I want to remove **B** from **C** to get **A**. **C** is a signal that is recorded by experiment, and so **A** and **B** are unknown. However, I can assume statistical properties of **A** and **B**. This problem has given me much headache because normally in the literature, **B** is treated as the desired signal to be extracted, and **A** is removed as unwanted "noise". I want to do this the opposite way around. This is the reason why I have tried (perhaps unnecessarily) to recast the problem in a different way, and this is also the reason why I have turned to statistical analysis. I need some guidance on which type of numerical procedure to use, and perhaps a good reference with some example problems. I suspect that some form of linear prediction algorithm used in numerical statistics might be beneficial here. * * * Here is more detailed information that might be pertinent. Given only **C** , but not **A** and **B** , I would like to approximate **A** and **B** using a statistical method. I am searching for an algorithm or method that is reasonably well-known and documented (i.e. a tutorial, book or paper is available, and the method is known to work well.) Here is what I know about the datasets. I am certain that these statements can be written in a more precise manner.   1. Datasets **A** and **B** are uncorrelated.   2. **A** and **B** are not constant.   3. **A** has smaller variability than **B**   4. Both the real and imaginary parts of **A** are "smooth" curves, where the rate of change is gradual   5. The real and imaginary parts of **B** are "spiky" curves, and the rate of change is quicker than **A**   6. The histograms of the real and imaginary parts of **B** are not skewed. Multiplying **A** by **B** skews the histograms of **C**. Can the histograms of **C** be separated in some way?   7. **A** represents the frequency-domain attenuation kernel of a signal, hence the need for complex variables.   8. **B** represents the signal in the frequency domain. The signal **B** is modified by **A**. The signal **B** is a reflection from a waveform that has been passed through an attenuating medium.   9. **A** models the effects of the medium on the waveform.   10. Despite **A** and **B** being in the frequency domain, they are nothing more than just datasets of complex numbers. Alternately, the problem can be reformulated so that **A** has elements consisting only of real values, and **B** is complex-valued, with real and imaginary parts. Thus, the imaginary part of **A** is zero for all elements in the dataset. The same statements above hold for **A** and **B**. Does this make the problem easier? I wonder if PCA or Fourier-transform methods would be useful for this problem. I am not certain which class of statistical methods would be useful for this type of problem.


title: Analysis of participant's results to see who is the most different

text: What I'm trying to do is analyse results for an experiment I have conducted. An example set of results is               Case A   B   C   D      1   2   3   1   8      2  11  11  14   5      3   2   3   4   3      4  12   7   8   7      5   6   6   9   4      A, B, C, and D are different participants. I would like to identify which participant has the most difference in their results compared to everyone else.


title: Choice of an evaluation metric for a graph clustering algorithm

text: I have instances for which the only thing I know is 70% of the distance matrix. I know some of these points form groups of correlated points (each point of a group is "close" to every point of the same group). I want to find these clusters of "close" points. I chose a distance threshold. Each pair of points closer than this distance is linked by an edge. I get then a graph to which I apply a graph clustering algorithm. In order to select the 'best' parameters, I want to choose a quality metric. My objective is to regroup the points that form very interconnected groups. I first tried to compare the silhouettes of the clusters for different parameters, but I feel uncomfortable that 30% of the pairwise distances are ignored, and I am not sure this metric is appropriate with respect to my objective. Anyone has a more appropriate measure?


title: Is it possible to compute a covariance matrix with unequal sample sizes?

text: I'm not sure if this question is correct, but is there a way to construct a covariance matrix for two vectors that have different lengths? If so, how? And would it have a size of $(m+n) \times (m+n)$ (assuming the two vectors are of length $m$ and $n$)?


title: Introducing dummy variables to least squares estimation

text: I have a simple algorithm for identifying similar products based on product characteristics. When those characteristics are size, weight and price (ie continuous variables), a simple least squares strategy is effective. For a given product, pick as similar those products with the lowest sum of squared differences in product characteristics. My issue now is I'd like to introduce into consideration characteristics that can only be represented as dummy variables. An example might be "is a plasma TV". What I have now is a dummy variable, appropriately named because I appear to be too dense to figure out how to incorporate it into the algorithm. I realize a graph based or nearest neighbor approach would be a natural fit, but I'd rather keep things simple: LSE can be accomplished easily with some php, nearest neighbor would mean I'd have to pull in R or somesuch. Is what I'm asking possible?


title: Interaction Terms and Logit Models

text: I have read the previous discussions on interaction effects and main effects, and I have a question on the subject. I am running a destination choice model (multinomial logit), and I have one continuous variable (distance in meters) and several dummy variables accounting for given behavior (main activity at destination) which I am hypothesizing as having an interaction effect. In this particular model, the effect of Distance is as expected significant by itself, and also are significant the interactions between Distance and the dummy variables, however, when I try to include the dummy variables alone to get their main effects, the model cannot be estimated even after including only one dummy variable. I illustrate this below:   1. Distance+Distance*Dummy2+Distance*Dummy3 = Significant main and interaction effects (s.e. of interaction not estimated)    2. Distance*Dummy2+Distance*Dummy3 = Significant interaction effects    3. Distance+Dummy1+Dummy2 = Model cant be estimated    4. Distance+Dummy1+Dummy2+Distance*Dummy2+Distance*Dummy3 = Model cant be estimated  So my question would be:   1\. if it is possible to account only for the interaction terms without the main effects given what I explained before?


title: Uncorrelatedness implies independence under some condition?

text: How would I go about proving the following statement: Given two random variables, $X$ and $Y$, if $E(h(X)g(Y)) = E(h(X))E(g(Y))$ for all bounded and continuous $h$ and $g$, then $X$ and $Y$ are independent. I have no idea how to begin the proof. Please provide any pertinent references.


title: Taleb and the Black Swan

text: Taleb's book "The Black Swan" was a New York Times best seller when it came out several years ago. The book is now in its second edition. After meeting with statisticians at a JSM (an annual statistical conference), Taleb toned down his criticism of statistics somewhat. But the thrust of the book is that statistics is not very useful because it relies on the normal distribution and very rare events: "Black Swans" don't have normal distributions. Do you think this is valid criticism? Is Taleb missing some important aspects of statistical modeling? Can rare events be predicted at least in the sense that probabilities of occurrences can be estimated?


title: Propensity score matching after multiple imputation

text: I refer to this paper: Hayes JR, Groner JI. _"Using multiple imputation and propensity scores to test the effect of car seats and seat belt usage on injury severity from trauma registry data."_ J Pediatr Surg. 2008 May;43(5):924-7. In this study, multiple imputation was performed to obtain 15 complete datasets. Propensity scores were then computed for each dataset. Then, for each observational unit, a record was chosen randomly from one of the completed 15 datasets (including the related propensity score) thereby creating a single final dataset for which was then analysed by propensity score matching. My questions are: Is this valid way to perform propensity score matching following multiple imputation ? Are there alternative ways to do it ? For context: In my new project, I aim to compare the effects of 2 treatment methods using propensity score matching. There is missing data and I intend to use the `MICE` package in R to impute missing values, then `twang` to do the propensity score matching, and then `lme4` to analyse the matched data. **Update1:** I have found this paper which takes a different approach: Mitra, Robin and Reiter, Jerome P. (2011) Propensity score matching with missing covariates via iterated, sequential multiple imputation [Working Paper] In this paper the authors compute propensity scores on all the imputed datasets and then pool them by averaging, which is in the spirit of multiple imputation using Rubin's rule's for a point estimate - but is it really applicable for a propensity score ? It would be really nice if anyone on CV could provide an answer with commentary on these 2 different approaches, and/or any others....


title: Proving a non-stopping time

text: Let me begin by first confirming that this is indeed the correct place to post this (other ideas I had were math.SE). That said, * * * Let $X_n$ be a Markov chain on the state space $\mathcal S$ and for $ y \in \mathcal S$ let $T_y = \min\\{ n \ge 1 : X_n =y\\}$ be the first return time to $y$. Let $W_y = T_y - 1$ be the time just before the first return to $y$   * Explain why $W_y$ is not a stopping time    * Show that the Strong Markov Property does not apply to $X_n$ at random time $W_y$.  **My Work** When showing that $W_y$ is not a stopping time, is it sufficient to write $$W_y = \bigcap_{i = 1}^{n-1} \\{X_i \ne y\\} \cap X_n = y$$ and claim that since $X_n$ does not belong to the set $\\{X_0, X_1, \dots, X_{n-1}\\}$, we have that $W_y$ is not a stopping time? * * * Then, for showing that the Strong Markov Property does not apply, can I write $$\mathbf{P}(X_n = y \mid W_y = n-1, X_{n-1} = i, X_{n-2} = x_{n-2}, \dots, X_0 = y) = 1 \ne p(i, y)$$ where $p(i,y)$ is the one step transition probability from $i$ to $y$?


title: Approximate multidimensional sequences

text: I'm trying to find some method for approximating sequences. I found Keogh's SAX symbolic approximation but it works for real valued time series and I have a database of about 3000 sequences, each one with a different number of elements, each element being of the form (a,r,b,c) where a, b and c are categorial variables and r is a real variable. It is an heterogeneous multidimensional time series and I'm trying first of all to find state of the art techniques for treating them. Please any hint on how to approach this? Alessandra


title: CDF in a discrete simulation (Jerry banks book)

text: Consider the experiment of tossing a single die. Let $X$ be number of spots on up face of die after toss. Then range space of $X$ is $R_x = \\{1,2,3,4,5,6\\}$. The discrete probability distribution for this experiment is: \begin{array}{ccccccc} x_i & 1 & 2 & 3 & 4 & 5 & 6 \\\ p(x_i) & 1/21 & 2/21 & 3/21 & 4/21 & 5/21 &6/21 \end{array} Now the author has given a table for the above experiment with its cdf (Cumulative distribution function), but I don't understand how it is produced. The table is \begin{array}{ccccccc} x & (-\infty,1) & [1,2) & [2,3) & [3,4) & [4,5) & [5,6) & [6, \infty)\\\ F(x) & 0 & 1/21 & 3/21 & 6/21 & 10/21 & 15/21 & 21/21 \end{array} How is the $F(x)$ value calculated ?


title: How to deal with neural networks when the data is unbalanced (but unlimited)?

text: If you assume that you have a plentiful/unlimited supply of training data, and you are using a neural network for classification, how do you deal with a situation in which a very high proportion of the data belongs to one class? For example, if you are doing binary classification and 99% of data comes from one class and only 1% the other. When you train the network using backpropagation, it normally uses the overall average error. If you keep the data with the same skewness as it comes in, will that cause it not to learn the rarer class properly? Another option is to discard data (presumably at random) from the more common class, until you have an equal proportion of data from each class, e.g. it is 50/50 in my binary example. But if you balance the data like this, then the network is being trained on data that is not representative of real life. The assumption that training data is plentiful (or unlimited) means that balancing the data need not reduce the size of the training set. If extreme skewness is inadvisable, is there a threshold at which this occurs, e.g. is it ok for it to be 30/70 but not 1/99 ?


title: PCA on Binary Data

text: I having binary data set (yes/no), so can I apply PCA on that. Is it mathematically correct to do that. In my opinion Binary variable can only be subjected to logical operations, so how it can be legal to apply PCA on the binary varaibles . Previously Question similar questions have been asked How to perform principal components analysis (or factor analysis) on binary data (using SPSS)? but I could not clear my doubt.


title: Does it matter than the difference between two parameters on the logit scale doesn't map to their difference on probabilty scale?

text: I am trying to understand some hypothesis testing results I am getting from a logistic regression analysis. Say I have 2 estimates of the log-odds of an event happening (the betas from my GLM). On the logit scale these these are               log.odd.trt1 <- -0.803410     log.odd.trt2 <-  0.772911      On the probability scale these are               prob.trt.1 <- 0.3092966     prob.trt.2 <- 0.6841503      The difference between the two log odds is -1.576321               diff.log.odds <- -0.803410 - 0.772911      If you convert this difference to the probability scale you get 0.1713171               inv.logit(diff.log.odds)      In contrast, The difference between the probabilities is -0.3748537               prob.trt.1 - prob.trt.2      -1.576321 is not -0.3748537. These seems reasonable given how the logit transformation works. However, my concern is whether that these means it is valid to conduct hypothesis tests that ask "is the difference between two log odds 0." This is in contrast to the question "is the log odds different from 0". Some context: I am using logistic regression to compare 2 treatments implemented over 6 years. I am interested in determining for each year whether there are differences between the 2 treatments. I am using the R package multcomp with a contrast matrix to carry out these pairwise multiple comparisons and also conduct a test for linear trend within each treatment. I am uncertain, however, if I am doing this correctly because the hypothesis that I am testing it essentially "is the difference between to logits zero" and the difference between two logits does not map directly to the difference between two probabilities. My inverse logit function is `inv.logit <- function(x) {1/(1+exp(-x))}`


title: How do you establish complete versus partial mediation in a simple mediational model?

text: I recently received this question from a student: > In a simple mediation model, if I have found the indirect effect (ab) to be > significant and the direct effect (c') to be small and insignificant, does > that mean I have full mediation or partial mediation? ![enter image description here](http://i.stack.imgur.com/KBESX.jpg)


title: How to compare the variance of two conditions with different means?

text: I have a series of vectors, which are measurements from one sample at two time points. First time point:   1.`[5,3,2,4,4,3,6,5]` (mean=4.000)   2.`[5,6,3,3,4,3,4,5]` (mean=4.125)   3.`[6,3,4,2,5,3,5,7]` (mean=4.375)   ... etc. grand mean=4.17 Second time point:   1.`[1,2,1,2,1,3,4,5]` (mean=2.375)   2.`[2,2,3,1,1,3,3,5]` (mean=2.500)   3.`[1,3,1,2,2,3,5,4]` (mean=2.625)   .... etc. grand mean=2.5 I want to see if the variance for each measurement/vector is significantly different between the two time points. However, the second measurement has a lower mean, which can therefore drive overall variance. How do you compare the variance of two conditions when the means differ?


title: Error bars in interaction plot for ANOVA

text: I notice whenever I see an interaction plot for say a simple two factor ANOVA there are no error bars present, just the points for the estimated means. Is it ever appropriate to display error bars in an interaction plot for an ANOVA? If yes when would you want to do this and how would they be calculated? If no why not?


title: Why is this normal probability plot graph skewed right?

text: First, I apologize for the bad image, but it's the best I can do with a mouse. Anyways, here it is: ![enter image description here](http://i.stack.imgur.com/D3WwL.png) So the data is skewed right, but the normal probability plot bends up and over what would be the approximate linear equation. According to different websites1, 2, a normal probability plot with data skewed right goes under the approximate line of best fit, and my graph looks more like it's skewed left. What am I misunderstanding? 1: itl.nist.gov/div898/handbook/eda/section3/normprp4.htm   2: www.basic.northwestern.edu/statguidefiles/probplots.html#Data%20Skewed%20to%20Right


How can I calculate the confidence interval of a mean in a non-normally distributed sample? I understand bootstrap methods are commonly used here, but I am open to other options. While I am looking for a non-parametric option, if someone can convince me that a parametric solution is valid that would be fine. The sample size is > 400. If anyone could give a sample in R it would be much appreciated.


title: Sample-size calculations for Benjamini-Hochberg, Westfall-Young, Holm-Bonferroni methods

text: I'm currently setting up an experiment where we might want to do multiple comparisons (i.e. comparing several treatments with a control at the same time). It is pretty straightforward to calculate the needed sample-size using the Bonferroni correction. However, the Bonferroni method is pretty conservative, and I'm worried that we're wasting time or resources getting more samples than we actually need. Are there ways to calculate needed sample-size for other correction methods, such as Benjamini-Hochberg, Holm-Bonferroni, Westfall-Young correction? Or, in your experience, are you likely to see any significant decrease (more than 5%) in sample-size using any of these other methods at all? The test in question is a simple comparison of treatment effects on a categorical outcome variable, with expected mean at 50%.


title: Why does bootstrapping the residuals from a mixed effects model yield anti-conservative confidence intervals?

text: I typically deal with data where multiple individuals are each measured multiple times in each of 2 or more conditions. I have recently been playing with mixed effects modelling to evaluate evidence for differences between conditions, modelling `individual` as a random effect. To visualize uncertainty regarding the predictions from such modelling, I have been using bootstrapping, where on each iteration of the bootstrap both individuals and observations-within-conditions-within-individuals are sampled with replacement and a new mixed effect model is computed from which predictions are obtained. This works fine for data that assumes gaussian error, but when the data are binomial, the bootstrapping can take a very long time because each iteration must compute a relatively compute-intensive binomial mixed effects model. A thought I had was that I could possibly use the residuals from the original model then use these residuals instead of the raw data in the bootstrapping, which would permit me to compute a gaussian mixed effect model on each iteration of the bootstrap. Adding the original predictions from the binomial model of the raw data to the bootstrapped predictions from residuals yields a 95% CI for the original predictions. However, I recently coded a simple evaluation of this approach, modelling no difference between two conditions and computing the proportion of times a 95% confidence interval failed to include zero, and I found that the above residuals-based bootstrapping procedure yields rather strongly anti- conservative intervals (they exclude zero more than 5% of the time). Furthermore, I then coded (same link as previous) a similar evaluation of this approach as applied to data that is originally gaussian, and it obtained similarly (though not as extreme) anti-conservative CIs. Any idea why this might be?


title: k-means quick cluster analysis error in SPSS

text: I am doing a k-means cluster analysis for 6 component factors obtained by PCA. For centers and seeds, I have carried out a hierarchical cluster analysis. When I run the k-means cluster analysis, i get a warning message: > Text: aggregate, The file referenced on the FILE subcommand does not have > the same number of splitfile variables as does the working file. My cluster center file (aggregate obtained by the hierarchical analysis), contains 6 variables and 5 clusters (exactly same as input to k-means analysis). I have also changed the first column title of aggregate file to cluster. Anything else I need to do to resolve this issue?


title: Measures of variable importance in random forests

text: I've been playing around with random forests for regression and am having difficulty working out exactly what the two measures of importance mean, and how they should be interpreted. The `importance()` function gives two values for each variable: `%IncMSE` and `IncNodePurity`. Is there simple interpretations for these 2 values? For `IncNodePurity` in particular, is this simply the amount the RSS increase following the removal of that variable? I greatly appreciate any enlightment :)


title: What if I know for sure that my target variable is not normally but Beta distributed?

text: As stated in the question, what sort of modeling technique would be most appropriate? Thanks


title: Conditional restricted Boltzmann machines on a time series dataset

text: **Preamble of the problem** I am currently trying to apply Conditional Restricted Boltzmann Machines on a time series dataset problem, in particular, the dataset constitutes of `10 day` stock market recordings, and in each day, `50 stock value` recordings are logged that are 5 minute apart from each other. Therefore, to my presumption, each day can be considered as a sequence `{x0,x1,...,x50}` and each recording in a day is a sample `x`. **My presumptions - please correct me if I'm wrong**   1. Each day of recordings is a sequence   2. The purpose of C-RBM is to map a sequence - a day of stock recordings - into a set of features that can be used by any regression algorithm as input   3. A different `RBM` is trained on each sequence **Question** How do we feed the sample recordings of a sequence to an `RBM` ? Do we concatenate all the features of a sequence into one large vector and train the `RBM` on that? Or do we use each sequence as a miniature dataset on which the `RBM` trains, but then the time element is not emphasized here? The question basically boils down to how, from an abstract point of view, can we use `C-RBM` to extract hidden features from sequences? Your suggestions will be highly appreciated. Thank you


title: Threshold and root mean square error for data-clustering

text: I am working on new incremetnal data-clustering method. I start with finding sum of entire data series (all attributes), then find sum of individual attributes, compute p-value by dividing one series sum with all attribute sum for both series. Then compute error by using sum of series, p-value and (1-p) value, using square root in denomintor. I would like to know 1) what is the significance of finding square root in error formula? 2) How to compute threshold which decides members of cluster, after computing closeness-factor? This algorithm begins with two data series S1 and S2. To generalize this algorithm later, we will call data series as DSi where i varies from 1 to n. j is the point in any DSi. The next step is find T, T(j)=S1(j)+S2(j) The probability of the outcome S1 can then be calculated as: Then compute p=sum of S1(j) / sum of T(j) Error is computed as = [p*T(j) - Si(j)] / sqrt (T(j) _p_ (1-p)) weight of every attribute = sqrt of (T(j)) closeness-factor = sum of [error*error(j) * weight (j)] / sum weight(j)


title: Classification with mislabeled data

text: I have a classification dataset where roughly 20% (maybe more) of the labels are incorrect. There is no way to know which labels are incorrect and no way to eliminate them in the future when further data is collected. A method I see in dealing with this is to train an ensemble of classifiers and then take only the training data which matches the ensemble majority vote. Are there any other algorithms/methods that are more resilient to data that is not labeled 100% correctly? Can we even treat this data as supervised learning? Is there anyway to trust the trained model or performance metrics such as accuracy and F1 score? Thank you for the help.


title: Fitting data in a generalized Pareto distribution and parameter estimation

text: I have log(return) data as time series, how I can fit this data in a Generalized Pareto distribution and estimate the parameters of this distribution, any kind of resource pointer with clear code would be good for me, If it can be done using open-source like R or Python that would be great.


title: Fusing Aggregate (Macro) and Micro Data?

text: I have a microdata recipient survey and a macro (aggregate) donor. How can I fuse the (binary categorical) data? The Statistical Matching techniques/software I'm familiar with are micro-to-micro. Furthermore, a web search revealed nothing (but perhaps my terms were wrong). One (possibly over-simplified/naive) approach is:   1. Build a probability distribution of the donor's response in question, broken out into (harmonized) common variables like age, income, etc...   2. A respondent in the recipient that most closely matches donor profile X (e.g.: Nearest Neighbor) is assigned a "yes" response with probability X (based on #1).  Any pointers, comments on the above, references and software (R is a good one) are greatly appreciated.


title: Sum of products in an expected value

text: _A box contains $n$ balls numbered from 1 to $n$. Suppose you take a ball at a time, putting it back on the box, until you pick a ball twice. How many balls are you expected to take from the box?_ Let $X$ be the r. v. of interest. Its support is every natural number from $2$ to $n+1$. If $k$ is one such number, to get the first repetition on the $k$-th pick, the $k-1$ previous ones must be distinct and $k$-th equal to some of those, therefore, $$\begin{align} \forall k\in\mathbb N\cap[2,\,n\\!+\\!1]\qquad \mathbb P(X=k) &= \frac nn\frac {n-1}n\ldots\frac {n-(k-2)}n \;\cdot\; \frac{k-1}n =\\\ &=\frac{(n-1)!}{(n-k+1)!}\frac{k-1}{n^{k-1}} =\\\ &= \frac{k-1}n\,\prod_{l=0}^{k-2}\left(1-\frac{l}n\right)\quad, \end{align}$$ so that $$\begin{align} E(X) &= \sum_{k=2}^{n+1}\,k\cdot\frac{k-1}n\,\prod_{l=0}^{k-2}\left(1-\frac{l}n\right) =\\\ &=\frac1n\sum_{k=2}^{n+1}\,k(k-1)\prod_{l=0}^{k-2}\left(1-\frac{l}n\right)\quad. \end{align}$$ However, the book answer is $$E(X) = 2 + \sum_{k=1}^{n-1}\prod_{l=1}^k\left(1-\frac ln\right)\quad.$$ What am I doing wrong? Or, if the answers are actually the same, how to show that?


title: Changepoints in R

text: I have the following dataset:               results <- data.frame(Date = c("A", "B", "C", "D", "E", "F", "G", "H","I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S"),                       P1 = c(0.43, 0.45, 0.57, 0.15, 0.5, 0.33, 0.26, 0.81, 0.43, 0.48, 0.14, 0.26,-0.21, 0.27, 0.37, 0.33, 0.68, 0.15, 0.44))      I want to know, if there are statistically significant changes of my observations. So I thought to use changepoint analysis. First I used the bcp-Package, with the following code:               c <- bcp(results$P1)     plot(c)       However, there are no changepoints according to this plot. Then I used the "changepoint" package and the following code:               var=cpt.var(results$P1, method="PELT")     plot(var)      Here I get three possible changepoints, but not were I supposed them to be (for example not in M, but in N). Can anybody explain me why? Or is there another way to show, if the values changed significantly from one observation to another?


title: Comparing two means and need to know if they are statistically different

text: I have the net revenue from 81 companies in one year and I'm comparing to the net revenue of the same 81 companies 3 years later. How do I compare (simple subtraction of the means or t-test?), and how do I know if it's statistically significant?