title: Validation of clustering results

text: I have a data which contains several columns which I later reduced using a PCA algorithms to two different components. I then applied the k-means algorithms to the data.   Now, how can I verify that my data clustered well into each group? Or how do I determine misclassification rate? For instance, using R, if I check the cluster vector say k$cluster against the labels of the data I had previously before clustering can I just draw a confusion matrix from that and assume that 1 in the clustered vector is equivalent to 1 in my labels?               col3    col2     Col1   lables                                                123     2.32      2.50    0                124    2.81      3.10     1          125    2.72      3.09     2          126    2.92      3.03     3          127    2.32      2.95     4           Please note this is a hypothetical data; my data is way bigger than this.


title: Robust parameter estimation for Exponentially modified Gaussian distribution

text: I'd like to test how well my data can be modeled by an Exponentially modified Gaussian distribution (Wikipedia) or Normal-exponential-gamma (NEG) Distribution. However, the parameter estimation (which involves Skewness) is not very robust when there are outliers in the data set. I've good experiences with Median based parameter estimation on this data set (see the related question Estimating parameters of a normal distribution: median instead of mean? ). Do you know a Median-based or similarly robust parameter estimation for the ExGaussian / NEG distributions? I've already tried trimming my data set, the results were clearly better afterwards. Yet, this obviously introduces a bias, that I would need to correct somehow.


title: Sum of multiple loaded dice

text: I have $n$ dice with $m$ sides. The $i^{th}$ dice will show value $0 \leq x_i \leq m-1$ with probability $0 \leq D_i(x_i) \leq 1$. What is the probability that the sum of the dice equals $\alpha$ Is there some approximation for $P(\alpha)$


title: Exact logistic regression, R, and elrm: data formatting and output

text: I use SPSS, but am forced to use R for exact logistic regression. So I'm brand new to R (and hating it so far) and also new to logistic regression. I've read the original elrm paper and looked at examples of its use. However, I can't find information on the questions below (after the data description). The fit of two models of cognitive processing was compared for each subject in each of 3 conditions. My binary dependent variable is whether the difference in model fits was significant or not (my "Success" variable below). I have three experimental Conditions: 0, 1, and 2. 0 is my reference group. My question is: is there an overall effect of Condition? If so, which conditions differ? The specific alternative hypothesis is that the proportion/probability of "success" should be greater in conditions 0 and 1 than in condition 2. My data look like ![original data](http://i.stack.imgur.com/OAFtP.gif) ...and so on. SPSS actually creates the dummy variables for you on the fly but they are easy enough to create explicitly. **Question 1:** I have read that to use elrm you have to enter the data such that the response variable represents success/number of trials. And as far as I can tell elrm doesn't create dummy variables automatically. I've seen examples of tables representing this data structure, but can't find any step- by-step examples of getting raw data into that format, espescially given a one-variable 3-levels situation. Is there an example out there that I'm missing? If not, is this what the data should look like? ![reformated data](http://i.stack.imgur.com/aVxOL.gif) I'm not sure how I'd enter the dummy variables into the formula...just as separate variables? **Question 2:** I can see how I can get the tests of the coefficients of the dummy variables. But I can't figure out how to get a test of the overall effect of the independent variable. I need to evaluate the overall effect of Condition before looking at individual conditions. Is there a way to get that out of elrm? (I found an example of this done for the aod package which runs regular logistic regression but not exact logistic regression.) **Question 3:** I can't find a description of what the p-value for individual coeffeicients represents in elrm. Is this is for the Wald test?


title: How to do multiple regression with limited experience and (hopefully) excel?

text: I am doing a study of how legal need relates to a number of predictors. Outcome Variable: Legal Need (Yes or No)   Possible Predictors: Age, Gender, Race, Ethnicity, Language, Clinic, Insurance Status.   Predictors based on bivariate analysis: Age, Language, Clinic, Insurance Status. All ORs were around 1.6 and were significant, typical 95CI between 1.2-2.0.   Number of subjects: 8803 I suspect that the one variable that is going to fall out of a multiple regression analysis would be Clinic, as the Asthma Clinic (vs the well-visit clinic) has in general, older children who don't speak English and are more on Medicaid (the other predictors from bivariate analysis). Knowing all of this, I still can't figure out how to run a multiple regression in Excel. I keep getting Pearson's coefficients like 0.06, which it then says is significant. Either way I think I need to go beyond this, maybe a stepwise forward or backward analysis? It's difficult for me to grasp when the outcome is a 0 (no) or 1 (yes), this doesn't really result in anything linear when the percent "yes's" are typically in the 15-30% range. Am I missing something? Besides years of stats training?


title: Analysing change between two variables measured at 3 points

text: I have two variables measured concurrently at 3 time points, let's call them wealth (`W1`, `W2`, `W3`) and anxiety (`A1`, `A2`, `A3`). Suppose that wealth and anxiety are uncorrelated. Now, I'm interested whether _change_ in anxiety is related to _change_ in wealth. What would be the recommended way of testing this? Naively I would just create 4 new variables `W2-W1`, `W3-W2`, `A2-A1` and `A3-A2` and look at correlation coefficients between `W2-W1` vs. `A2-A1` and between `W3-W2` vs. `A3-A2`. Somehow this approach seems silly, since I want a single model that would look at the overall change. Should I be looking at mixed/multilevel models? Specifically which model? Thanks!


title: Criticism, please: simplifying state space model

text: As a little experiment, I am extending a nice, interpretable AR/MA relationship between a security $r$ that is variably influenced by the previous $k$ time points over another security $f$. These relationships ("lags") are represented by $l$. I want to pump up the "realism" by making both securities move with stochastic volatility: $ \begin{bmatrix} r_t \\ f_t \end{bmatrix} = \begin{bmatrix} l_{r,1} & ... & l_{r,k} \\ l_{f,1} & ... & l_{f,k} \end{bmatrix} \Theta_t + \Sigma_t u $ $u \sim N(0,1)$ $ \log(\Sigma_t) = \begin{bmatrix} \log(\sigma_{r,t}) \\ \log(\sigma_{f,t}) \end{bmatrix} = \begin{bmatrix} 1 & \phi_r \\ 1 & \phi_t \end{bmatrix} \begin{bmatrix} \omega_r & \omega_f \\ \log(\sigma_{r,t-1}) & \log(\sigma_{r,t-1}) \end{bmatrix} + \begin{bmatrix} \eta_{r,t} \\ \eta_{f,t}) \end{bmatrix} $ $\eta_{i,t} \sim N(0,\sigma_{\eta_i})$ $ \Theta_t = \begin{bmatrix} f_t \\ f_{t-1} \\ ... \\ f_{t-k} \end{bmatrix} = I \begin{bmatrix} f_{t-1} \\ f_{t-2} \\ ... \\ f_{t-k-1} \end{bmatrix} + \begin{bmatrix} \sigma_{f,t-1} \\ \sigma_{f,t-2} \\ ... \\ \sigma_{f,t-k-1} \end{bmatrix} $   1. Does this look right to better trained eyes, before I whip out a prototype? Especially, suggestions on parameters that might to be worth consolidating (i.e., params not adding much explanatory "value")?    2. Any suggestions for how to approach the "lag" selection more rigorously (e.g., model selection procedure)?    3. Does anyone see shortcuts for consolidating the two transition equations ($\log(\Sigma_t)$ and $\Theta_t$) into the typical "state space representation," i.e., the $\Theta$, $F$, $G$, $V$, $W$ referred to in the literature. For some special cases the equations can be simplified substantially: e.g., traditional MA log volatility models, a log transformation of the master equation simplifies out that pesky $u$ stochastic term.


title: How to understand complex social relations in huge populational data?

text: I have 20 years of data from observing dolphins. When a group is seen, it gets an unique number identifying it, and all identified (marked) dolphins were also registered. So I had a table like               Group Dolphin     1     1     1     10     1     14     2     10     2     23      I got about 20,000 groups and 500 dolphins. I made a matrix with groups in rows and dolphins in columns, with 1 where a dolphin is in a group, 0 where it's not. So I've ran a multivariate analysis (prcomp in R), and it gave me two very distinct blocks. ![PCA](http://i.stack.imgur.com/CjplT.png) I know that the vertical axis (PC2) is related with the size of the group, the small groups are below, in the bottom of the triangles. But the two blocks are separated in the horizontal axis (PC1), and I didn't figured yet what it's about. I suppose the dolphin population might be separated into two "clans" or something like that, but couldn't discover yet. I colored the points according to group size, water level, place of sight, but all patterns appear equally on both sides. The block on the right has about 2400 points. My questions are: such a separation must mean something, right? What other tools can I use to understand what is going on? Thanks very much!


title: How to equate questionnaires with different response scales

text: I'm trying to help out a not-for-profit organization and made an error that I'm trying to fix. They've been doing a non scientific survey of their Board of Directors for a number of years regarding board engagement. They asked me to help them convert the survey format from an unfriendly Excel file to something web based, handle the results and do some analysis for them. In past years, they've used what I'm going to call a weighted scale (Very good = 5, Good = 3, Fair = 2, Poor = -1). This year, we mistakenly included an extra level by allowing a 5, 4, 3, 2, or 1 response. (That is, the previous questionnaires had 4 options, but this one has 5.) I'm trying to figure out how to conform the scales to be roughly equivalent. I had been thinking that I could let 5 = 5, 4 = 3, 3 = 2.5, 2 = 2 and 1 = -1. However, I realize that this doesn't work and will skew everything lower. Can anyone recommend a good way to handle this?


title: Comparison of two populations by analyzing the dependence of a single dependent variable on several independent ones

text: I am looking at two populations on which I have measured 4 independent variables ($X_1$, $X_2$, $X_3$, $X_4$) and one dependent variable ($Y$) along with each measurement or observation. I suspect that the two populations are different in the sense that the second population's $Y$ shows a dependence on an additional non-considered factor (I expect $Y$ to be systematically lower, in fact). Considering that:   * the dependent variables are highly correlated;   * the dependence of $Y$ on the $X$'s is unknown and arguably non-parametric;   * the number of observables/measurements is small (~100 for the first group, ~30 for the second). I feel that multi linear regression would assume some kind of underlying structure and would bias my analysis, so I would like to avoid it. I tried to use PCA on both and compare the PC's: nothing significant comes out in the sense of dimensionality. Could you suggest a test that would check if the two populations are really different in respect to $Y$ while keeping at a minimum the number of assumptions on its functional dependence? Thank you very much, P.


title: How to model zero inflated, over dispersed poisson time series?

text: I am trying to model weekly disease counts in 25 different regions within 1 country over a ten year period as influenced by temperature. The data is zero inflated and over dispersed. I am most familiar with Stata but I don't think that there is any option amongst the `gee`, `xtmixed`, `xtmepoisson` etc. commands that allows me to account for the zero inflation and over dispersion issues as well as the autocorrelation. I log transformed the incidence data and used a SARIMA model but the residuals are not quite normal. I think that there are versions of the ARIMA model for integer data like disease counts but I can't find a program for it. I was also thinking that I could create a hierarchical model with random intercepts for each region and random effects of temperature in each region, while also accounting for the regular seasonal disease cycle. I believe that I could model this in R using a package like glmm.admb but due to my limited statistical and R knowledge I am not entirely sure how to do use it. I am mainly confused about accounting for the autocorrelation and seasonal cycle part of the data using a program like this. Any advice on how to best do this?


title: Computing a Gaussian posterior from a Gaussian prior and likelihood function in R

text: I'm new to both R and Bayesian statistics, and I have a problem where I have a normally distributed prior that elicits a mean and standard deviation. The introduced likelihood function is also normally distributed with a mean and standard deviation that can be drawn from a sample. Now I understand that a posterior is formed that will also be normally distributed. I have been asked to generate this in R, but I cannot find an example of R code where the posterior is formed from a normally distributed prior and likelihood function. Could someone please point me in the right direction. Thanks.


title: PCA: 91% of explained variance on one principal component

text: I am new to PCA and wanted to do a bit of experimentation on my data set just to see what it looked like (using R). I am not able to give access to the data here since it is confidential. However, if there is some other kind of statistic/visualization you would like to see that would help you answer my questions please let me know and I will provide it. I found the following information about the explained variance:               Component Prop.Var     1         0.911804348     2         0.033618098     3         0.020827269     4         0.011772988     5         0.006611746     6         0.005372772     7         0.004464788     8         0.003436401     9         0.002091589      This raises the following questions:   1. Am I justified in removing the other 8 principal components?    2. How do I interpret 91% of explained variance on one component?   3. If I only kept one component what would be the best way to visualize the data? Below is how the graph of the first two principal components looks. The spread of the data like this is not surprising given how little of the variance is on the second component. ![Principal Components 1 and 2](http://i.stack.imgur.com/g9wXP.png) As I mentioned, I am new to PCA so I really do not know if there is even any useful information to be found from this kind of dimensional reduction. Any insight would be appreciated.


title: Why Baum-Welch algorithm is an instantiation of EM algorithm?

text: $\newcommand{\E}{\mathrm{E}}$ I don't understand why Baum-Welch algorithm is an instantiation of EM algorithm. Indeed, why computing $\alpha_t(i)$ and $\beta_t(i)$ corresponds to Expectation step. Expectation step corresponds to compute expectation over the latent variable of log-likelihood of observed variable given the parameter: $\E_{Z|X,\theta^{t-1}}[L(X|\theta)]$. I don't understand the link.


title: Analyze results from a/b test using welch's t test

text: I'm performing a web based a/b test where there is a control and one treatment. The results are not as simple as "converted" or "didn't convert." A user can "convert" anywhere between 0-10 times. I have all of the data from both the control and the treatment (how many times each user converted from both the control and the treatment). The sample sizes aren't the same size (about 1900 vs 2100) and the variances are different (7.12 vs 6.02). The mean of the treatment is about 11% higher than the control. The goal of the experiment is to find out if the treatment can increase converts per user. The numbers are showing an 11% increase in conversions per user. To find out if the result is statistically significant I've been trying to use a Welch's t test. When I use the equation for a Welch's t test found on this wikipedia page I get the following results:               t-score: 2.26     degrees of freedom: 4025.82      On the wikipedia page it says I can use a `t-distribution` to test the null hypothesis (my null hypothesis is that the means of the control and treatment are the same). But I'm not sure how to go about using a `t-distribution`. I'm guessing it has something to do with using a t distribution table. Questions:   * Is using a Welch's t test a good approach for this situation?      * If not what other method would you suggest?   * With my results from the Welch's t test, how do I use a `t-distribution` to determine if I have a significant result and with how much confidence?


title: Joint expectation of normal distribution

text: I have a vector $v=(v_1,...,v_J)$, which is jointly normal. I then need to take the expectation of this vector. This should give me something like this: 1)$$E(v)=(E(v_1),...,E(v_J))$$ However, I want to write up the integral for this expectation, but I'm unsure how to do this. My guess is something like this: $$E(v)=\int v \phi(v)dv$$ How does this relate to the expectation I have in 1)? $\phi(.)$ is the joint pdf, but to compute something like $(E(v_1),...,E(v_J))=$, I guess I need the marginal distribution: $$E(v_1,...,v_J)=(\int v_1 \phi(v_1),\int v_2 \phi(v_2),...,\int v_J \phi(v_J))$$ Thanks


title: Very small coefficients

text: Suppose that you are running a linear regression model with outcome $Y$ and explanatory variables $X_1$ and $X_2$. You also run a linear regression model with a quadratic transformation of $X_1$: $$E[Y] = \beta_0 + \beta_{1}X_{1}+\beta_{2}X_{1}^{2} + \beta_{3}X_{2}$$ If $|\beta_2| \approx 0$ can we drop it from our model?


title: Another manner for handling Tobit - related regression in R

text: Which method would you suggest for estimation an univariate regression ($y = a + bx + e$) which comprise a lot of $0$ values in the dependent variable (actually kind of censored data). My experience hints at the use of the **Tobit model**. But There is a problem connected with the normal distribution and homoscedasticity of the error term. In my case I'm suspecting a _violation_ of these basic assumptions for using Tobit. In particular the behaviour of latent $y$ under the zero is not anticipating. Do you know some method (literature) for handling this kind of problems?


title: t.test in R between two normal distributions obtained by using the dnorm function

text: I would like to calculate Student's t test of two normal distribution data. However, I do not have the data at all, but only mean and standard deviation of each population. So, how can I simulate those data in R and calculate the t.test among both? These are the values (mean and SD): Population 1: 6,62 +- 0.52 years Population 2: 6.31 +- 0.49 years Thanks in advance Mario


title: How to initiate bias node in a restricted Boltzman machine

text: I am new to Neural Networks and trying to implement RBM. I am stuck on initializing the visible layer's bias value. Is it supposed to initialize to some random number or there is some probabilistic distribution by which it should be initialized. In python I am seeing:               > "vbias = theano.shared(value =numpy.zeros(n_visible,     >                               dtype = theano.config.floatX), name='vbias') "      but don't know what is happening here in this line.


title: Predictor is significant in logistic regression, but not in Mann-Whitney

text: I have run a logistic regression using several predictor variables (call them p1 p2 p3 p4) to predict a binary dependent variable (call it y). P1 is a significant predictor of y in the regression. However, a Mann-Whitney test of p1 shows no significant difference by category of y. Is it normal to have a significant logistic regression predictor that doesn't significantly differ by groups, or is there potentially something screwy with the data?


title: Comparing Gini coefficient of 2 distributions

text: I wonder if it is fair is to compare the gini coefficient of two discrete distributions with different number of elements. If not, how can I adjust the coefficients for a fair comparison. In particular, I'm interested in the case when the sum of values in both distributions are equal. For example, minutes played by players of two teams in a basketball game if they have different number of players (10 players in one team, 12 players in the other team) but for both teams sums plays 48 minutes x 5. Other example, different partitioning of some surface in different regions, how to compare inequality for different partitioning if they have different number of subareas.


title: Model instability in data mining. When it is big enough to discredit a model and how to measure it?

text: Let's say I have two models. One has cumulative lift on test data 4.322578, second 2.84488. The only advantage of the second over the first consists in the quality of having the cumulative lift curve nearly identical for both train and validate sets. For the first model the curve of cumulative lift looks as you can see below. ![enter image description here](http://i.stack.imgur.com/GtD7X.png) Sorry for non-english elements on the plot, but that's the only SAS Enterprise Miner version I have access to. Mean standard error for the model from the plot - the first one - for train, validate and test sets are: 0.038389, 0.04999 and 0.055314. For the second one: 0.070393, 0.071889 and 0.07727. First model seems to be obviously better in terms of accuracy. But at the same time it is less stable, more sensitive to data changes. My question is whether those stability issues are here big enough to discredit that model in favor of the second one? How to measure instability and at which point it renders a model unusable? If it is important, the models in question are neural networks.


title: R plot set scale

text: I did some research and didn't find what I was looking for: I have a data set x = 10,100,1000 and y1, y2, y3 respectively, and I want to plot those points. But what is happening is that in the x axis the values that are appearing are 10, 200, 400, 600, 800 and 1000. How can I set the scale so that it only appears the values I specified above? What I'm doing is:               y_values <- c(9, 9.6, 9.7, 9.8)     x_values <- c(10,100,1000)     plot(x_values, y_values)      Thank you.


title: Asymptotic convergence involving EDF

text: Please help me proving this: Suppose that $Y_1,\ldots,Y_n$ are i.i.d. nonnegative RV's with CDF $F$ and $E(Y_i)=\mu<\infty$. Let $y_1,\ldots,y_n$ be a realization from which an EDF $\hat{F}$ is defined: $$ \hat{F}(y)=\frac{1}{n}\sum_{i=1}^nI(y_i\leq y). $$ **Want to show** : $$ \text{plim}_{n\to\infty}\left(n^{1/2}\int_0^\infty y(\hat{F}(y)-F(y))d(\hat{F}-F)(y)\right)=0. $$ Please state any further necessary assumptions. * * * EDIT responding to GUNG'S comment: This result is a small step in a bigger argument in Russell Davidson 2009's paper on the Gini index. The author must have considered it obvious because he mentions it without proof. I have some understanding of asymptotic theory such as some laws of large numbers, some variants of the CLT, various modes of convergence, Slutsky Theorem, measure theory in general and probability theory in particular, and some stochastic calculus. My problem is I've been learning asymptotic theory on my own and I haven't done a lot of exercises in this. So while I want to understand the author's reasoning, I don't know where to start with this one.


title: Significance test for PMCC

text: 149 students each sat two tests. My null hypothesis is that the product-moment correlation coefficient is 1, because each test should be an accurate indicator of current attainment. The actual PMCC for the data is 0.912152856. Please could you show (and explain) how to test if this is a significant result?


title: GBM, it's overfitting/multicollinearity problem and parameter setting up

text: I recently came across a predicting problem (0-1 outcome, with more than 80 variables), I decided to use GBM (Gradient Boosting Machine by Friedman)to handle this job. I let the GBM use only 70% of the dataset, and within this training set, 60% are used at each iteration. During this process I found a big problem of GBM: there was one predictor (AAA) that was so important (nearly causal) that this single predictor can predict the outcome with precision=40% at recall=40% (0-1 unevenly distributed, with 0: 80%, and 1: 20%), and this single predictor model was better than the GBM with 80 predictors including it (AAA). I wonder why this could happen, I would like to hear you thought on it. From a 'linear regression' perspective, I suspect multicollinearity, that GBM splits many other spurious predictors that correlated to AAA, gradually missing the predictive power of AAA. Additionally, I found GBM cut each of my predictor in a very 'overfitting' way (predictor sometime positive, sometime negative affect the fitted value). To tackle this problem, my strategy was to use a parameter of GBM called 'var.monotone' (it takes three values, 1 for positive correlation, -1 for negative, 0 for arbitrary). I selected 22 out of the original 80 predictors based on there association with the outcome, and used var.monotone to force their relationship with the outcome. the model was slightly better. I would like to know if force each predictors relationship with the outcome is a good strategy to mitigate overfitting, what's the better way from your perspective? More importantly, do you think a blackbox with 80 predictors really work? how about a logistic regression, or a L1 regularized one, with only 5 highly important (or causal) predictors? Thanks, Hicham Hu


title: Statistical measures for data validation

text: I have to validate the data in one database by comparing it to data from a validated application system. The solution I came up with was (at least in theory) creating a statistical method that would define the number of random samples needed to check on a table of the database, in order to say, that with a 95% chance the number of errors is representative for that table. Then, with the probability of the errors I could consider if the data is valid or not. The binomial distribution doesn't help, since I don't have any knowledge about the probability at which an error could be found. My statistical knowledge is pretty rusty and I can't figure it out right now. Just to clarify again: Given N the number of entries in a table I want to know what number n for random samples I have to chose in order to be able to say after testing, that the probability of the error is p, at a 95% accurateness. I would be very helpful for every helping idea I could get. Any references to literature or sites covering this would also be greatly appreciated. Thanks in advance!


title: Two ways of using bootstrap to estimate the confidence interval of coefficients in regression

text: I am applying a linear model to my data: $$ y_{i}=\beta_{0}+\beta_{1}x_{i}+\epsilon_{i}, \quad\epsilon_{i} \sim N(0,\sigma^{2}). $$ I would like to estimate the confidence interval (CI) of the coefficients ($\beta_{0}$, $\beta_{1}$) using bootstrap method. There are two ways that I can apply the bootstrap method:   1. Sample paired response-predictor: Randomly resample pairs of $y_{i}-x_{i}$, and apply linear regression to each run. After $m$ runs, we obtain a collection of estimated coefficients ${\hat{\beta_{j}}}, j=1,...m$. Finally, compute the quantile of ${\hat{\beta_{j}}}$.   2. Sample error: First apply linear regression on the original observed data, from this model we obtain $\hat{\beta_{o}}$ and the error $\epsilon_{i}$. Afterwards, randomly resample the error $\epsilon^{*}_{i}$ and compute the new data with $\hat{\beta_{o}}$ and $y^{*}_{i}=\hat{\beta_{o}}x_{i}+\epsilon^{*}_{i}$. Apply once again linear regression. After $m$ runs, we obtain a collection of estimated coefficeints ${\hat{\beta_{j}}}, j=1,...,m$. Finally, compute the quantile of ${\hat{\beta_{j}}}$.  My questions are:   * **How are these two methods different?**   * **Under which assumption are these two methods giving the same result?**


title: Use of R for biostatistics publications

text: Presumably this has been discussed already, but I wanted to get some sense of how prevalent is the use of R for biostatistics publications. On PubMed I did find a number of research papers that were conducted using R (as mentioned in the Methods or Statistical Analysis sections of the documents). However, the majority were done using SAS/SPSS and also an appreciable number using Stata. There has been a lot of talk about R in recent years and I am looking for some data that can corroborate the fact that R has been rising in popularity in the commercial area. In my experience, however, I have been seeing some resistance for firms to switch to R from SAS, even though there is so much more that you can do with R compared to SAS. There is also the concern that FDA "prefers results in SAS" ... although I have not seen anything officially mentioned by FDA to this effect. I understand that SAS has been there for decades, it is what many R&D departments are most comfortable working with, etc., etc., ... and am looking for real-world examples of who are instead using R for professional research and/or publications.


title: Can I use AIC value for comparing logit and probit model where for each model the number of covariates are equal?

text: Hi i am trying to use AIC value for comparing logit and probit model where in each model the data and the number of covariates are same (say, covariate= 3 for each model) Does AIC value increases while increasing sample size? Which one better to use as a tool for comparing logit and probit, AIC or MSPE?


title: When y= 1, logit is infinity. How can you regress that? Yet somehow, that's logistic regression

text: Logistic regression essentially means taking the logit of your response proportions, p, and then doing standard regression. Consider the case where one p is 1. logit(1) = log(1/(1-1)) = infinity. How can you do regression with infinity? That is, if one of your observed proportions p is 1, then you are trying to find the line that minimises the sum of squared differences from a set of points that include infinity. How does this not drive the regression line to always have infinite slope?


title: Randomly select one data entry for the morning and one for the evening for each day in R

text: I have a data set of recorded elephant locations for the span of 6 years (elephantdata). Some of the entries are duplicates of the same day. I want to create a subset of my 'elephantdata' and create a new matrix that only has one data point before noon and one data point after noon for each day. I would like that data point to be selected randomly from the data points offered. For example:               elephantdata     Row     Date       Time    Breeding_H   E1  E2  E3      1     11/01/06    08:38    1   0   0   0      2     11/01/06    18:00    1   0   0   0      3     12/01/06    05:38    1   0   0   0      4     02/02/06    08:58    1   0   0   0      5     02/02/06    16:30    1   1   1   1      6     02/02/06    17:15    1   1   1   1      7     03/02/06    08:45    1   0   0   0      For Breeding_H I want to find that if the 'Date' is a duplicate AND two entries are in the morning (before 12:00) or in the evening (after noon) I want R to randomly chose one of the entries. So Row 6 and 7 are both after noon and on the same day, so I would like to only have one of those data entries. Any suggestions? Would really appreciate any help you've got. I've looked up the 'unique' and 'duplicated' function but I don't know how to write the actual code to make this happen.


title: Self-study question on Confidence Interval (true/false)

text: Good evening everyeone, I am currently working on a set of self-study questions which relate to true/false answers. I am currently faced with a question where the answer claims for it to be false, but I find it true. Not sure if there is an error in the answer. Appreciate some help and advice please. > A confidence interval is an estimate for which there is a specified degree > of certainty that the sample statistic will be in the interval. My take is that the above **statement is true** as the confidence interval tells us with a specific degree of confidence how likely a value will fall within the range. Appreciate some advice please. ![enter image description here](http://i.stack.imgur.com/mzKSx.jpg)


title: ABC. How can it avoid the likelihood function?

text: According to the Wikipedia article, we have the scenario shown below, but how can ABC generate simulation datasets from samples of $\theta$ without knowing or evaluating the likelihood function? For example, how would this method work for a normal $N(\mu, \sigma)$ distribution? Assuming that the parameter to estimate is $\mu$ and that $\sigma$ is known, how would it generate the simulations without evaluating $N(\mu, \sigma)$? ![enter image description here](http://i.stack.imgur.com/pkpG0.png)


title: Cluster analysis

text: I am trying to cluster cells (1×1km) over a specific area. Each cell is composed of various habitats defined by a code. (Each habitat consists of 3 parameters, so a habitat code looks like e.g. 1-3-15. There are around 100 different habitats, i.e. max number of combinations of these 3 parameters). I now try to define clusters according to:   * Number of habitats per cell   * 1st largest habitat covering the cell   * 2nd largest habitat covering the cell and   * Rarest habitat in the cell (I developed a global rarity index for each habitat, and one habitat having the lowest index is chosen per cell). So the last 3 in the list refer to a specific habitat per cell, and the first is more a quantitative measure of how many are in each cell. Based on this they'd like to perform a cluster analysis, but I am wondering whether there is a mix up between categorical and quantitative data. Could anyone please give some advise on what to method to use, e.g. in R? Thanks.


title: What is the difference between dprobit and probit in stata?

text: My boss gave me this do file but I have never worked with dprobit. I also need to interpret it, this is the result:                xi: dprobit apoyado $X, robust nolog     i.fuerza          _Ifuerza_1-9        (_Ifuerza_1 for fue~a==Comando General  omitted)     note: _Ifuerza_8 != 0 predicts success perfectly           _Ifuerza_8 dropped and 1 obs not used          note: _Ifuerza_2 dropped because of collinearity     note: _Ifuerza_9 dropped because of collinearity          Probit regression, reporting marginal effects           Number of obs =    592                                                             Wald chi2(5)  =  33.04                                                             Prob > chi2   = 0.0000     Log pseudolikelihood = -387.38545                       Pseudo R2     = 0.0433          ------------------------------------------------------------------------------              |               Robust      apoyado |      dF/dx   Std. Err.      z    P>|z|     x-bar  [    95% C.I.   ]     ---------+--------------------------------------------------------------------     _Ifuer~3*|  -.1759987   .1586211    -1.10   0.273   .162162   -.48689  .134893     _Ifuer~4*|  -.3915721   .1314311    -2.56   0.011   .155405  -.649172 -.133972     _Ifuer~5*|  -.2404857   .1450899    -1.60   0.110   .533784  -.524857  .043885     _Ifuer~6*|   .0252801   .1673778     0.15   0.881   .087838  -.302774  .353335     _Ifuer~7*|   .0776492   .1834521     0.41   0.681   .038851   -.28191  .437209     ---------+--------------------------------------------------------------------       obs. P |   .5675676      pred. P |   .5723459  (at x-bar)     ------------------------------------------------------------------------------     (*) dF/dx is for discrete change of dummy variable from 0 to 1         z and P>|z| correspond to the test of the underlying coefficient being 0


title: Create a threshold for binary classification problem based on distribution of criterion

text: I created a criterion for my financial data-set for classify data to two classes for other processing (like `neural network binary classificatio`). After calculating this criterion (between 0 and 1) the distribution of output is like `exponential` (but between `0` and `1` - something like `beta distribution`). Now i want set a `threshold`to create two classes from these samples. I can set every threshold (`70% - (0.7)` was good in `classification accuracy`) but I want i want a logic for that. for example in `normal distribution` we can set `68–95–99.7 rule`. I want set this default value in my program because different samples are inserting in my system in the feature. Thanks.


title: Prior selection for Gaussian Processes (GP)

text: I am trying to select a prior for the covariance parameters of my Gaussian Process (GP) and have been running into numerical problems with my MCMC code. My model is the following: $$Y = D\beta + GP(0,\sigma^2R(\psi))$$ where $R$ is the following correlation matrix: $$R_{i,j}(\psi)=\exp(\psi(x_i-x_j)^2)$$ My priors for $\beta$ and $\sigma^2$ are the following: $$p(\beta)\propto 1$$ and $$\sigma^2\sim IG(a,b)$$ so what I am trying to figure out is what are possible appropriate choices of prior for $\psi$. So far I have been trying to place an inverse gamma prior on $\psi$ but that has not been going very well. Here is an update: Here is further problem formulation if it helps Likelihood: $$f(\beta,\sigma^2,\psi|x)\propto|\sigma^2R(\psi)|^{-1/2}\exp(-\frac{1}{2}(Y-D\beta)'R(\psi)^{-1}(Y-D\beta))$$ Prior: $$p(\beta,\sigma^2,\psi)=p(\beta|\sigma^2,\psi)p(\sigma^2|\psi)p(\psi)\propto\sigma^{-a-1}\exp(-b/\sigma^2)p(\psi)$$ Now, based on my posterior under this prior, I have been trying to code up a Metropolis within Gibbs sampler with the following full conditions: $$\beta|\sigma^2,\psi,Y\sim N(\hat\beta,D'(\sigma^2 R(\psi))^{-1}D)$$ $$\sigma^2|\beta,\psi,Y\sim IG(a,b)$$ $$\psi|\beta,\sigma^2,Y\propto |R(\psi)|^{-1/2}\exp(-\frac{1}{2\sigma^2}(Y-D\beta)'R(\psi)^{-1}(Y-D\beta))p(\psi)$$ So I sample $\beta$ and $\sigma^2$ within a gibbs step and then use a M-H step for sampling from $\psi$ but I am having horrible numerical problems when sampling from $\psi$ (basically my acceptance ratio is 0/0). Finally, I'll post some code to if that helps. Code:               x = sort(runif(5,0,1))     y = exp(-1.4*x)*cos(7*pi*x/2)          plot(x,y,type='l',ylim=c(-1,1))          tau = as.matrix(dist(x,upper=T,diag=T))          corrR = function(psi,tau){         R = exp(-tau^2*psi)         return(R)     }          psi.cond = function(psi,r,sig,D,beta,Y){         d = .01         e = .01         ans = det(r)^(-.5)*exp(-.5/sig*t(Y-D*beta)%*%solve(r)%*%(Y-D*beta))*psi^-(d+1)*exp(-e/psi)         #ans = -.5*log(det(r))-.5/sig*t(Y-D*beta)%*%solve(r)%*%(Y-D*beta)-(d+1)*log(psi)-e/psi         return(as.real(ans))     }          D = rep(1,length(x))     Y = as.matrix(y)          m = length(x)     a = .01     b = .01          B = 100     beta = rep(NA,B)     sigma = c(1,rep(NA,B-1))     psi = c(120,rep(NA,B-1))          for(i in 1:B){              R = corrR(psi[i],tau)         bhat = as.real(solve(t(D)%*%solve(R)%*%D)%*%t(D)%*%solve(R)%*%Y)              beta[i] = rnorm(1,bhat,t(D)%*%solve(sigma[i]*R)%*%D)         sigma[i+1] = 1/rgamma(1,(m+2*a)/2,(as.real(t(Y-D*beta[i])%*%solve(R)%*%(Y-D*beta[i]))+2*b)/2)                   log.xi = rnorm(1,log(psi[i]),.1)         xi = exp(log.xi)              u = runif(1)              R.xi = corrR(xi,tau)         R.psi = corrR(psi[i],tau)              temp = (psi.cond(xi,R.xi,sigma[i],D,beta[i],Y)*(1/psi[i]))/(psi.cond(psi[i],R.psi,sigma[i],D,beta[i],Y)*(1/xi))         alpha = min(1,temp)              if(u <= alpha){             psi[i+1] = xi         }else{             psi[i+1] = psi[i]         }          }


title: Whether to perform a 3 by 2 repeated measures ANOVA or three separate paired t-tests?

text: I have $2\times 3$ repeated measures design.   * 2 conditions(a and b)   * 3 intensities (low, medium and high). I am concerned only with the a to b differences at the corresponding intensity. For example I would like to compare:   * A low and b low   * A medium and b medium   * A high and b high. I am not concerned with comparing A low to A medium to A high as physiologically we know that they are different and it does not answer my research question. I don't care about how the 2 variables interact, but I may be wrong. My inclination is to say the paired T -tests are fine as **Should I perform 3 paired t-tests for the comparisons of interest or a $2\times3$ Repeated measures ANOVA?**


title: How can two positive dependent correlation coefficients differ significantly without differing significantly from zero?

text: I have 3 variables A, B and C, which are all markedly skewed. Below, I indicate correlations between two variables, e.g., A and B as $r_{ab}$. Using the following test with with $n-2$ degrees of freedom $$t= \frac{r\sqrt{n-2}}{1-r^2}$$   * $r_{ac}$ is positive but not significantly different from zero,    * $r_{bc}$ is positive but not significantly different from zero, and    * $r_{ab}$ is significantly different from zero. Then using Williams t test   * $r_{ac}$ and $r_{bc}$ are not significantly different from each other. How can this be possible? **Shouldn't the difference between two positive correlations always be non-significant when they are both not significantly different from zero?**


title: Connection between PDFs/PMFs and Bayes Theorem

text: **UPDATE** Original question was confused and poorly worded. I thought about it more and don't think I have a question any longer. After thinking a bit more I came up with: For a distribution, such as the normal for example, the probability that the variable falls in a certain range is:               P(a<=X<=b|X ~ N(mu, sigma^2))       Using Bayes' Theorem I can infer:                P(a<=X<=b|X ~ N(mu, sigma^2)) = P(a<=X<=b, X ~ N(mu, sigma^2)) / P(X ~ N(mu, sigma^2))      Interpreting this statement I could say: P(X ~ N(mu, sigma^2)) is the prior probability that the random variable X ~ N(mu, sigma^2) (based on whatever information we have), and P(a<=X<=b, X ~ N(mu, sigma^2)) is the prior probability that X ~ N(mu, sigma^2) **and** within the interval [a,b]. Then, if it is revealed to us that X is in fact distributed N(mu, sigma^2), or if we assume it to be so, the probability that it falls within the interval [a,b] is the ratio of the prior probabilities, as stated above. The above statement doesn't depend at all on what the particular pdf is that you're dealing with. I think I was hoping to see a connection between the two, but they seem to be unrelated ideas.


title: Random walk with bivariate normal distribution

text: Let $X$ be a random variable from $f(x; \theta)$, where $\theta =(\theta_1,\theta_2)$. I want to: take a sample from this distribution using Metropolis Hastings algorithm and update the parameters simultaneously. The proposal distribution is a bivariate normal distribution with mean the current value of the parameter $\theta$ and covariance matrix $\Sigma$. Steps: Given $\theta$ Step 1: Generate $\theta^{can}$ from the proposal distribution. Step 2: Take $\theta$ = $\theta^{can}$ with probability $\alpha(\theta,\theta^{can}) $ where $\alpha(\theta,\theta^{can})$ = $ min \left( \frac{f(\theta^{can})}{f(\theta)}∗ \frac{proposal(\theta)}{proposal(\theta^{can})}, 1\right)$ I know that $proposal(\theta^{can}) = \frac{1}{\sqrt{(2\pi)^k|\boldsymbol\Sigma|}} \exp\left(-\frac{1}{2}({\theta^{can}}-{\theta})^T{\boldsymbol\Sigma}^{-1}({\theta^{can}}-{\theta}) \right)$. I would like to ask if $proposal(\theta) = 0$. Thank you very much.


title: Data formatting for an ANCOVA in SPSS

text: was wondering if someone could help clarify an ANCOVA for me, or maybe post me in the right direction. I have a set of outcome measures that have been administered pre treatment and post treatment. I do not have a control group. My data is in a long format (i.e. each individual has data on 2 rows, one for pre treatment one for post treatment). I want to do an ANCOVA with baseline scores as a covariate. Time is my fixed factor. I've been told to have my data in the long format (by someone quite senior), but I'm wondering whether I need to transform it into wide format, and then use outcome score as my DV. I'm using SPSS for the analysis.


title: Sample size from graph

text: I have a screen shot and the question asks if the sample size is blocked out, what is the most likely sample size that was used to generate these results. How do I do this? The possible answers are 5, 10 , 20 and 25 and 25 is the correct answer. It gives the mean 8.08, median 7, sd=6.22 skew=0.83 at the top.


title: What type of statistics is used to describe this chart?

text: What kind of data analysis technique / formula would you use to describe which is the fastest OS in this chart, Windows or Linux? ![enter image description here](http://i.stack.imgur.com/uleOH.png) If Linux and Windows had equal speed on all tests, all bars would be 50% blue, 50% green. The data is duration in milliseconds. Each number (e.g. Linux Startup) represents the mean of 1000 individual duration tests on that scenario. The results of each test individually is covered in separate charts. This data represents an overview of the one set of data (Linux over the course of these 5 tests) and how it relates to the second set of data (Windows over the course of these 5 tests). The aim is to be able to say: Over the 5 tests we Windows is faster. (You can see this in the chart because there is less green, but this is to prove it statistically). Note: The chart below shows the same data in real value terms, showing that one test dominates the others. While this is an important result, it does not weigh on the current question as the two results in the 'Refreshing Desktop' test still bear the same relation to each other as those in the other tests. In the comments there has been some discussion on how to represent this data. The reason I represent it in stacked bars is because I want to show it _without_ the overshadowing effect of the one anomalous data. Is there a preferred method to show this relationship? ![enter image description here](http://i.stack.imgur.com/NqQNC.png)


title: Is there any decision tree version where each node selected by some measures over histograms other than single numbers?

text: What I mean is, constructing a decision tree where each node is not a single feature instead composition of multiple features, therefore the evaluation criteria should consider a metric over multiple numbers instead of single number. For the original case of decision tree we are looking for a feature's value comparing its ordinal relation < or > to other instances. But I am expecting to compare multiple value like $x_1,\cdots,x_n$ with new instance's values $x_1,\cdots,x_n$ at each node selection. It is like to have an vector of values selected for each decision-tree node.


title: Computing BIC for SUR model

text: Consider the following m regression equation system: $$r^i = X^i \beta^i + \epsilon^i \;\;\; \text{for} \;i=1,2,3,..,T$$ where $r^i$ is a $(T\times 1)$ vector of the T observations of the dependent variable, $X^i$ is a $(T\times k)$ matrix of independent variables, $\beta^i$ is a $(k\times1)$ vector of the regression coefficients and $\epsilon^i$ is the vector of errors for the $T$ observations of the $i^{th}$ regression If the above is an SUR model (seemingly unrelated regressions), does it make sense to compute the BIC (Bayesian Information Criterion) for the model? If so, how would that be computed? Should the dependent variables be stacked into one vector, in order to get only one regression equation? Any help would be appreciated! Thank you!!!


title: Is it true that adjusted-$R^2$ not a measure of fit? Why or why not?

text: Previously I've read that adjusted-$R^2$ is not a measure of fit. Recently, though, I wanted to substantiate that piece of knowledge by understanding the reason why but I couldn't find any substantive sources to back this up. The Wikipedia article on it states "while $R^2$ is a measure of fit, adjusted $R^2$ is instead a comparative measure of suitability of alternative nested sets of explanators" but does not provide a citation. http://en.wikipedia.org/wiki/Coefficient_of_determination


title: Differences in AUC calculation in R between pROC and AUC

text: I was comparing the performance of pROC and AUC libraries when performing auc() calculations on random data:               library("pROC");     library("AUC")          predictor <- rnorm(10000, 5);     outcome <- rnorm(10000) > 0;          print ("pROC:::auc() time & output")     system.time(x <- pROC:::auc(outcome, predictor))     print(x);          print ("AUC:::auc() time & output")     system.time(x <- AUC:::auc(AUC:::roc(predictor, factor(outcome))))     print(x);      AUC:::auc seemed to perform substantially faster, but what I found strange is that the compute different auc values for the same dataset:               > system.time(x <- pROC:::auc(outcome, predictor))        user  system elapsed         1.00    0.01    1.31           > print(x);     Area under the curve: 0.5058          > print ("AUC:::auc() time & output")     [1] "AUC:::auc() time & output"          > system.time(x <- AUC:::auc(AUC:::roc(predictor, factor(outcome))))        user  system elapsed         0.19    0.00    0.18           > print(x);         [1] 0.4942452      I thought the Auc() function was deterministic so they should produce the same number. Yet pROC produces 0.5058 and AUC produces 0.4942452 . Am I misusing either function? EDIT: FYI I tried making the number semi random and the functions now give identical results (bar rounding errors):               predictor <- runif(10000);     outcome <- as.integer((predictor + runif(10000)) > 0.5);


title: How to use the Pareto distribution in fitting survival curves?

text: I have a series of numbers, which are some survival probabilities that form a decaying curve. I would like to fit them with a "Pareto" distribution. I expect to have a smooth fitted curve, thus I can extrapolate it and have a rough prediction. My data are like `1, 0.987, 0.972, 0.965, ...`. How can I achieve my goal?


title: Statistical significance seen with Kaplan-Meier but not with Cox model

text: I am currently working on a retrospective analysis where I am looking at patients to see if immune sensitization prior to a procedure has any bearing on outcomes after the procedure. The immune sensitization is a categorical variable (Hi vs Lo). In regards to outcomes, one of the things I am looking is survival analysis (time to death). I ran a Kaplan-Meier analysis broken down by immune sensitization and there seems to be a statistically significant (p<.05) difference in regards to the time to death. My questions are these:   1. I ran a Cox regression analysis with the sensitization status as the only covariate and noted the p-value to be .35. Why is there such a big difference in the statistical significance and which one should I choose? (if it helps, from a visual inspection, only patients in the Hi category died)   2. When analyzing my patient characteristics, I noted that there was a gender difference in the two sensitization categories (p<.05). I ran a Cox Regression model with the sensitization and gender as the covariates and noted that p-value of 0.955 for the sensitization status and 0.5 for the gender. Am I correct in interpreting that neither covariate has any bearing on the time to death? * * * _As requested, here is some of my output_ Original Kaplan Meier analysis (PRA_Categ represents the sensitization status): ![Original Kaplam Meier Results \(PRA_Categ represents the sensitization status\)](http://i.stack.imgur.com/ba1aJ.png) The Corresponding COX regression with only one covariate (PRA_Categ): ![The corresponding COX regression with only one covariate \(the PRA_Categ\)](http://i.stack.imgur.com/4Eabs.png) COX regression with two covariates (PRA_Categ and Gender): ![COX regression with two covariates \(PRA_Categ and Gender\)](http://i.stack.imgur.com/R2Bht.png)


title: Controlling for variables in negative binomial regression

text: I would like to run a negative binomial regression, but also controlling for one of the variables. Similar to what you would do with a hierarchical regression. I am trying to predict violence (never, once, 3 times ...); and the literature says I should control for age, as aggression decreases with age. Would welcome your thoughts, thanks in advance. sharon


title: Derivation of the effect of unmodeled confounders on OLS estimates

text: I'd like to see how not adjusting for confounders affects estimated coefficients in linear regression. Here is what my lecture notes say:   1. Un-adjusted model: $E[Yi] = \beta_0 + \beta_1X_i$   2. Adjusted model: $\quad\ E[Yi] = \gamma_0 + \gamma_1X_i + \gamma_2W$, where $W$ is the confounder. Let $\hat{\beta_1}$ be the OLS estimate of $\beta_1$. Then under the adjusted model we have: $$ E[\hat{\beta_1}] = \gamma_1 + r_{XW}\sqrt{\frac{{\rm var}(W)}{{\rm var}(X)}}\gamma_2 $$ **How $E[\hat{\beta_1}]$ was calculated?**


title: $p$-values from two-sample $t$-tests of resampled vector clump at $p > 0.7$

text: I have created a vector of length 400 with 6 unique values of relatively equal proportions (1-6). I have drawn two random samples from the vector, each $n=26$, 10,000 times, performed $t$-tests, added the $p$-values to a list and output a histogram of the $p$-value frequency. The histogram of the $p$-values shows a cyclical clumping effect for $p$-values greater than $0.7$. I would like to know why this occurs. Here is my R-code:               set.seed(117)     cover1<-rep(1,100)     cover2<-rep(2,63)     cover3<-rep(3,59)     cover4<-rep(4,74)     cover5<-rep(5,56)     cover6<-rep(6,48)     landscape=c(cover1,cover2,cover3,cover4,cover5,cover6)          sample.size<-26     x<-10000     pval.list <-vector("list", length(x))     alpha <- 0.05     conf.level <- 1 - alpha     for(i in 1:x){            s1 <- sample(landscape, size = sample.size, replace = TRUE, prob = NULL)       s2 <- sample(landscape, size = sample.size, replace = TRUE, prob = NULL)       result.t.test <- t.test(s1,s2, conf.level=conf.level)       pval.list[i] <- result.t.test["p.value"]      }     s<-pval.list     pvalue<-as.numeric(s)     hist(pvalue, main = paste("Histogram of" , "p-value"), breaks= 100)      ![Histogram of p-values](http://i.stack.imgur.com/Hz4pH.jpg)


title: Does PCA preserve observations (rows) of the data?

text: Say I have a data matrix of size $N \times P$ where $N$ is the number of samples and $P$ is the number of features. Now, if I do principal component analysis, I get another data matrix of size $N \times K$ where $K$ was chosen according to some criteria. My question: if I pick a row (sample) from the $\text{PCA}$ matrix, does it still point to the _same_ sample as in the original data matrix? In my study the data on each row is from one subject, so I want to know if the correspondence still exists if I use $\text{PCA}$ for feature selection. (I think this is correct but better safe than sorry...!)


title: Odds that the jack of trump dealt 46% of time in pitch

text: In the game of Pitch what are the odds of the jack of trump being dealt in the 24 cards?


title: What are advantages of multi-point statistics compared to variogram?

text: **Variogram** is simple but effective way to investigated the spatial variation in the variable of interest in the field being studied. However it has been reported that variogram cannot be considered as the most satisfactory model especially where there are curvature _(e.g., permeable channels in reservoir engineering)_ in the medium being studied. That is different patterns could have more and less exact copy of the same variogram.   A suggestion for development an approach was **multi-point statistics** ( **geostatistics** ), in which for example, a **training image** , a 2D matrix template is being used for inference of spatial variation. It was shown in the literature that training images are useful to dictate the output the desired pattern honoring conditioning data and being based on stochastic approaches to satisfy statistical inferences.   * Questions: > **Q1:** What else about advantages and or disadvantages of multi-point > statistics ?   >  **Q2:** How to obtain training images?   >  __ **Q2.1:** How to generate training images using existing knowledge etc?   >  **Q3:** How about 3D and further dimensions?   >  **Q3:** What is the next stage after training image?


title: Time-varying linear regression

text: I've got two time-series, S (stress in a railway track) and T (temperature). The time-series are several months long. The relationship is linear, however, it can change subtly or 'jump' at some points, due to work done on the rails. Currently, I'm doing linear regression on a per-day basis, and then temperature-correcting S based on that days models. However, changes in the relationship do happen within each day, so this is not perfect. I guess I'm trying to figure out how to best fit multiple linear models between S/T over time. How would I go about doing that? I've just started out with R, though I'm also good with Python.


title: Chi-squared test for histogram data after doing an averaged shifted histogram

text: I have a data set of 903 continuous observations, that I graphically visualize with a histogram. The bin and width values could be optimized, but it is logical from the distribution that I have a Gaussian function. When I do the fitting, I use the frequency values of the data as the Y-values. For example, if the observations are ${(2,3,3,3,4,4,5)}$, and the user-defined bin width is 1.0, then the corresponding y-values would be $(1,3,3,3,2,2,1)$, assuming that the first bin limit will be assigned as the minimum value. I am not obtaining statistical significance with a Gaussian fitting, in terms of goodness of fit ($Q$) with the Chi-squared test ($\chi^2$). In other words, my null hypothesis is rejected making that a Gaussian model does not represent the experimental data. Now, I am doing the same test but using an averaged shifted histogram version of my x-values and the frequency for those averaged bins as the y-values for each observation, I now obtain good results in terms of goodness of fit. I need to clarify if it is valid to realize a Goodness of fit for an averaged shifted histogram, or if there is clear bias for data overfitting. Here it is a q-q plot of the data: ![QQ Plot](http://i.stack.imgur.com/gRHTp.jpg)


title: Restricted least squares: condition is that the parameter vector is greater or equal zero

text: I only found a example in which the constraind on $\beta_i, \ i=1,...k$ would make them sum up to one as well introduces a constraind such that all $\beta_i$ are all greater or equal zero. Since I dont need a condition for the summation of the $\beta_i$ I'm not quite sure how to deal with this problem. Since in this context a matrix $C$ of constraints is introduced such that $C\beta = c$, I'm not sure how to deal with $C\beta \geq 0$. Any suggestions are very welcome!


title: Is there any research on control of hidden stochastic processes?

text: **The short question:** Are there any studies of systems described by several hidden stochastic time dependent variables and observed variables that are given by known deterministic functions of hidden variables? **More Details:** I have a time dependent stochastic process in which $n$ real-valued variables depend on time: $x_i(t)$. These variables cannot be observed directly ($i.e.$ they are hidden variables). However, we can observe a deterministic (not stochastic) function of these hidden variables: $y(t) = y(x_{i}(t))$ (in my particular case $y(t) = \sum_{i=1}^n x_{i}(t)$). Moreover, we can observe another variable that is a function of hidden variables. This second function is also deterministic but it combines hidden variables corresponding to different times: $z(t) = x_{0}(t) + x_{1}(t-\Delta t) + x_{2}(t-2 \Delta t) + ... + x_{n}(t-n \Delta t)$. In addition to that we know that subsequent vectors representing the hidden variables are related in a stochastic way: $x_{i}(t) = \psi[x_{i}(t - \Delta t)]$. The exact form of $\psi$ is not known. We also always know $x_0(t)$ and we set the value of $x_{n}(t)$. The goal is to find the stochastic relation $\psi$ and to chose such values of $x_{n}(t)$ that $x_{t}=0$ (we choose $x_{n}$ now and what to keep $x_{0}=0$ in the future). We also want to avoid $y=0$. My question is how to approach this problem and if there are some results for such systems (or what are the name of such systems). The closest that I can think of are Markov hidden process but, as you can see, it is something different.


title: Can I trust the result of a one-way ANOVA with many treatments each with very few replicates

text: I am doing a one-way ANOVA of a response variable (Y) on a treatment factor (T) of 7 levels. However, for each treatment, I have only 2 observations (or replicates). The ANOVA result shows the the means across the treatment are highly significant (p<0.001). I have check the normal and constant-variance assumptions and seems that my data does not violate these assumptions. My question is: Can I trust the result of ANOVA, even if the sample size is extremely small (n=2*7=14)? The first question that comes to my mind is the power of my test. Since my ANOVA has obtained an effect size of 0.9. I conducted a power analysis (using Q*Power.3) to see how many samples are needed to obtain such effect size with significance level of 0.05 and power of 0.8. The result showed that I need at least 28 samples. Since I have much fewer samples than required, does it mean that my ANOVA results are not reliable?


title: GARCH-M(1,1) where ARMA(0,0) is "removed" in R

text: Which of the following is the correct code for fitting a GARCH-M(1,1) model where the ARMA(0,0) is "removed"? Or what is the correct code?   1. `modm11 = ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,1)), mean.model=list(archm=T, archpow=2, include.mean=F))`   2. `modm11 = ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,1)), mean.model=list(armaOrder=c(0,0), archm=T, archpow=2, include.mean=F))`   3. `modm11 = ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,1)), mean.model=list(armaOrder=c(0,0), archm=T, archpow=2))`


title: LaTeX output for R's summary.lm object - while displaying the information outside the table

text: This seemed to me to be basic, but I can't seem to find a solution online, so I wondered what I might be missing. I wish to include the output of an lm summary object inside an Sweave (.Rnw) document. I can either output the summary.lm as is, or use the xtable/Hmisc packages (through xtable or latex commands). Is there something like xtable that also gives the summary information which available from outside the table? ($R^2$, F statistics etc...?) Thanks.


title: Prewhitening time series: ARIMA-modelling versus polynomial trendelimination

text: I'm concerned with Box-Jenkins-models and especially the first step, the prewhitening to obtain meaningful crosscorrelations for identifying transfer functions and building regression models. I'm working with SAS and here one has to identify a suitable ARIMA-model for the input series to achive white noise and afterwards the response series is automatically filtered by the same model. Now I'm wondering about two things:   1. If I have two series A and B and I do not know which one is the input and which one the response series (and there can be feeback maybe), do I have to set at first series A as input and filter it with a suitable ARIMA-model (and the series B is filtered like that) and I can only have a look at positive lags (and lag 0) of the crosscorrelation function (CCF) that will show how series A affects B. And as a second step do I have to choose series B as input and do the same (set a suitable ARIMA-model for B, A is filtered with the same parameters, look at the positive lags (and lag 0) of the CCF to see how B affects A) ? It seems a little bit weird to me, but it seems to be right in a logical way. But then I ask myself what to do with two different lag 0 values for example? (In some cases I got significant results in one case and when I changed input and response like described it wasn't significant)   2. What if I decided that there is a deterministic trend in one of the series (let's say series A) and I have to detrend by polynomial subtraction? In this case in SAS the response series (let's say series B) isn't filtered in the same way. Is it right to take the detrended series from A and the original series B for computing the CCF? I would be very glad if someone can help me. Greetings, Daniela


title: What is a good probability model of the distribution of frequencies as represented by western sheet music?

text: I'm attempting to answer a general question about a method of musical comparision which is generalizable beyond western music but uses western music scales as a nearly sufficient subspace for dimensional reduction. What is a good probability model of the distribution of frequencies as represented by western sheet music for SATB and instruments?


title: How to decide the parent distribution of a sample?

text: 1) Suppose, I have 25 sample values. I have no idea as to which distribution, the observation comes from. How will I proceed? 2) With regard to the same question: Suppose, though I have no idea, I pretend that it comes from normal distribution and estimate its mean and variance, do a P-P or an Anderson-Darling and conclude that it comes from normal density. Precisely, does it essentially mean that the data follows normal distribution. Aren't I just trying to fit my data to a convenient distribution? In short, given a continuous sample, how do I decide its original parent distribution??


title: How do I know which class my random forest model is predicting?

text: I have a random forests model with which I am trying to predict species presence or absence.   This is my code:               #read in dataframe containing observations of species presence/absence & predictor variables     mydata <- read.csv('mydata.csv')          #fit random forests model     fitmodelA <- randomForest(SPECIESA ~ var1 + var2 + var3 + var4 + var5 +var6 + var7 + var8 +      var9 + var10, data=mydata, mytry=3, ntrees=500, replace=T, importance=T, keep.forest=T)          #predict to new data     predictmodelA <- predict(fitmodelA, newdata, type="prob")          # save as raster image     writeRaster(predictmodelA,"predictSPECIESA.tif")      Apparently I should get back a matrix that has the probability for both classes, i.e., in two columns. Do I understand correctly that it is also possible to add the "index" argument to predict only one class? With my code the way it is, my output raster produced one layer with probabilities 0 to 1, but no other attributes – **what class is this predicting?** I am more interested that my map show predictions of presence rather than absence. Probably a simple solution to this...? Thanks!


title: Using panel data with different variables: SUR or OLS separatedly?

text: I have a merged database from two cross-sectional surveys for distinct years ($t_1$ and $t_2$). If all questions were exactly the same, I could just stack the data horizontally and run a normal panel data OLS. The problem is that some variables exist only for $t_1$, while other variables exist only for $t_2$. So my question is: should I just do the same as above and disregard these variables, should I run cross-sectional OLS on each separate databases using all variables available for that year (and then there is a way to "aggregate" the results?) or could I run a SUR model where for each year I use all available variables and constrain the coefficients to be the same?


title: ARIMAX feature selection

text: I am implementing the ARIMAX model and need to implement feature selection. I have >100 features and a lot of data, so I need a method that isn't too computationally expensive. I tried a wrapper, backward feature selection algorithm but that was too expensive. Is there a less computationally expensive approach towards feature selection? Perhaps like a filter feature selection method.


title: Confidence Interval for squared term evaluated at the mean in Stata

text: When running an OLS regression with a squared term, $$ y = a + b_1(X) + b_2(X^2) + e $$ I know that the partial effect of $X$ is $b_1 + 2b_2(\bar X)$ to get the overall effect of $X$ on $Y$ evaluated at the mean. Using               reg y x c.x#c.x     margins, dydx(*) atmeans      does this in Stata and it also includes a 95% CI for the overall effect. My question now is: how does `margins` calculate the standard error and the 95% CI of the overall effect of $X$ at the mean?


title: Data processing before applying SVM

text: I am working on classification of audio files. It is a binary classification and I plan to use SVM. I have used SVM before for face matching and other image analysis and retrieval stuff. I have extracted the required feature vectors from the audio files, i.e., the training and test dataset and reduced their dimensionality by using Principal Component Analysis. I would like to know whether there are any more steps which are necessary before applying SVM classification and prediction? Should the test and training dataset obtained after applying PCA be normalized or centered? Would the results be different (better / worse) after applying the normalization / centerization? Or are there any more methods that can be used to pre-process data before SVM is applied on it?


title: Does the central limit theorem hold for the prediction error over different samples?

text: Given an (infinite) data population from which you repeatedly draw samples of a fixed size. On each sample you learn a classifier which you then evaluate by computing the prediction error on a large independent test set. The prediction error is defined as the average over all instances in the test set of the zero-one loss function $\mathcal{L}(y,\hat{y})$. (The zero-one loss function $\mathcal{L}(y,\hat{y})=0$ if the predicted label $\hat{y}$ of an instance equals the true label of the instance y, and is 1 otherwise.) My question is whether the central limit theorem holds in this situation, so that the distribution of the prediction error over the samples approximates a normal distribution, given a sufficient number of samples? I know that the theorem holds for a sufficiently large number of _independent_ random variables. So I am in doubt because the individual predictions for one test set originate from the same classifier, and the classifiers vary over the samples.


title: Which statistical test should I use for my experiment on aggressive interactions in killifish?

text: I am doing a project on sexual selection (male-male competition) in the turquoise killifish Nothobranchius furzeri. There are two morphs of male in the population from which my fish are obtained from- one has a red tail and the other has a yellow tail. My null hypothesis is: Tail colour is not related to dominance/competitive ability. I will be putting one yellow-tailed male and one red-tailed male in an arena and recording the number or aggressive interactions that take place within 5 minutes and the winner of each. I have 8 red males and 8 yellow males all of similar colour intensity. I originally thought I had about 30 of each and I was going to rank for size and pair them up (i.e. largest red with largest male etc.)which would control for size (the larger fish are more dominant). However with a sample size of only 8 this would not be sufficient to get significant results. How could I redesign the experiment and which statistical tests could I use? I don’t have access to any more fish.


title: Generating random samples from a density function

text: How does a computer algorithm set up to take as input an arbirary bivariate probability density function, generate pairs of numbers from that distribution? I have found a routine called simcontour that is part of LearnBayes in R that performs that operation.


title: Pseudo R squared formula for GLMs

text: I found a formula for pseudo $R^2$ in the book Extending the Linear Model with R, Julian J. Faraway (p. 59). $$1-\frac{\text{ResidualDeviance}}{\text{NullDeviance}}$$. **Is this a common formula for pseudo $R^2$ for GLMs?**


title: What is the difference between ITT and ATE?

text: I am having trouble understanding the different estimators that can be used in an impact evaluation. I know that the intention-to-treat (ITT) estimator compares differences between eligible individuals without the program, and eligible individuals with the program, regardless of compliance. However, I thought the average treatment effect (ATE) also measured the same thing. However, it seems that the ATE takes into consideration the compliance. Therefore, it compares outcomes between those eligible and taking up treatment with those who are not eligible. Is this correct?


title: Where can I read about the justification for the use of parametric probability distributions?

text: I would like to find a reference, preferably free on the internet, where I can read about the theoretical or practical justification for the use of parametric / analytic probability distributions. By parametric distributions I mean the named ones like Normal, Weibull, etc.


title: Probabilities in case-controlled studies

text: I have a nested-case control study that I have been using for analysis. At the end of my work I have deduced a set of variables that I use later to to classify new cases. One example of a simple classifier I am using is a naive Bayes, which will output simply a probability. So here is my question: Could I make my probabilities reflect the real world? In my specific example, the condition that I am testing for has a prevalence of 33% in my study, but a it has a population prevalence of only 10%. Bayes factors have been suggested to me as a way to achieve this, however I am little unsure how to set up the problem. As an example I have seen a Bayes factor as a logit between the true vs. study prevalence of the outcome. The classifier however was a logistic regression, and in that case the Bayes factor was just added to the linear predictors. I think the example there was very specific, and perhaps an inappropriate method for probabilities of a naive Bayes. Instead what I did was add the logit Bayes factor to the logged probabilities, but I am also not convinced this is right either. I also think a simpler solution would be to use Bayes theorem directly, but there I am not sure how to represented my study vs.population prevalences. The method below isn't quite right, but gets at what I want:                       p_final = classier_posterior*(population_prev)/(study_prev)      I should contextualize that I use the probabilities to establish a threshold for classification down stream.


title: What does "AR(p) filtered series" mean?

text: I guess this means that omitting some variables in a certain interval, say, $(x_1, x_2, x_3, x_4, x_5) \to (x_1, x_5)$ in AR(4) model. Is it right? Or does this means eliminating autocorrelations by adjusting autocorrealated variables by a certain method? Thanks in advance!


title: Transform odds ration (OR) in the context of multivariate meta-analysis in R (metafor package)

text: I would like to conduct a multivariate meta-analysis (multiple treatment arm meta analysis) comparing the effect of different drugs. My outcome measure is discrete and describes the number of occurrences of a specific side effect during the period of observation. The studies that i would like to include report the odds ratio (OR) and their 95%-CI of this side effect for every drug relative to a “baseline drug”. The problem is that studies defined different drugs to be the baseline.   * Is it valid to just transform my OR to bring them to a common baseline drug?    * Is it valid to transform the 95%-CI of the OR in the same way?   * There is typically no 95%-CI reported for the baseline drug. Thus   when I transform this baseline drug, I have no 95%-CI. Is there a way to derive a 95%-CI in this case? Here is some dummy data:               structure(list(study = structure(c(2L, 2L, 2L, 2L, 2L, 3L, 3L,      3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L), .Label = c("Liu", "Mark", "Smith"     ), class = "factor"), drug = structure(c(1L, 2L, 3L, 4L, 5L,      1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L), .Label = c("A", "B",      "C", "D", "E"), class = "factor"), OR = c(1, 1.5, 1.7, 1.8, 2.5,      2.8, 1.1, 1, 2.3, 1.2, 1.8, 1.2, 2.5, 1, 1.8), ci_lb = structure(c(13L,      5L, 4L, 9L, 11L, 12L, 2L, 13L, 10L, 1L, 6L, 3L, 8L, 13L, 7L), .Label = c("0.5556",      "0.6225", "0.7619", "1.2628", "1.365", "1.4212", "1.5758", "1.6208",      "1.7048", "1.7435", "1.9497", "2.3531", "NA"), class = "factor"),          ci_ub = structure(c(13L, 2L, 7L, 5L, 10L, 11L, 1L, 13L, 9L,          4L, 8L, 3L, 12L, 13L, 6L), .Label = c("1.5775", "1.635",          "1.6381", "1.8444", "1.8952", "2.0242", "2.1372", "2.1788",          "2.8565", "3.0503", "3.2469", "3.3792", "NA"), class = "factor")), .Names = c("study",      "drug", "OR", "ci_lb", "ci_ub"), row.names = c(NA, -15L), class = "data.frame")          > my_data        study drug  OR  ci_lb  ci_ub     1   Mark    A 1.0     NA     NA     2   Mark    B 1.5  1.365  1.635     3   Mark    C 1.7 1.2628 2.1372     4   Mark    D 1.8 1.7048 1.8952     5   Mark    E 2.5 1.9497 3.0503     6  Smith    A 2.8 2.3531 3.2469     7  Smith    B 1.1 0.6225 1.5775     8  Smith    C 1.0     NA     NA     9  Smith    D 2.3 1.7435 2.8565     10 Smith    E 1.2 0.5556 1.8444     11   Liu    A 1.8 1.4212 2.1788     12   Liu    B 1.2 0.7619 1.6381     13   Liu    C 2.5 1.6208 3.3792     14   Liu    D 1.0     NA     NA     15   Liu    E 1.8 1.5758 2.0242               my_data_1 <- my_data[my_data$study=="Mark",]         my_data_2 <- my_data[my_data$study=="Smith",]     my_data_3 <- my_data[my_data$study=="Liu",]          my_data_2$OR <- my_data_2$OR * 1/my_data_2$OR[my_data_2$drug=="A"]     my_data_3$OR <- my_data_3$OR * 1/my_data_3$OR[my_data_3$drug=="A"]     my_data_2$ci_lb <- my_data_2$ci_lb * 1/my_data_2$OR[my_data_2$drug=="A"]     my_data_3$ci_lb <- my_data_3$ci_lb * 1/my_data_3$OR[my_data_3$drug=="A"]     my_data_2$ci_ub <- my_data_2$ci_ub * 1/my_data_2$OR[my_data_2$drug=="A"]     my_data_3$ci_ub <- my_data_3$ci_ub * 1/my_data_3$OR[my_data_3$drug=="A"]          my_data_new <- rbind(my_data_1, my_data_2, my_data_3)     my_data_new     my_data_new$OR <- round(my_data_new$OR, 2)     my_data_new$ci_lb <- round(my_data_new$ci_lb, 2)     my_data_new$ci_ub <- round(my_data_new$ci_ub, 2)          > my_data_new        study drug   OR ci_lb ci_ub     1   Mark    A 1.00    NA    NA     2   Mark    B 1.50  1.36  1.64     3   Mark    C 1.70  1.26  2.14     4   Mark    D 1.80  1.70  1.90     5   Mark    E 2.50  1.95  3.05     6  Smith    A 1.00  2.35  3.25     7  Smith    B 0.39  0.62  1.58     8  Smith    C 0.36    NA    NA     9  Smith    D 0.82  1.74  2.86     10 Smith    E 0.43  0.56  1.84     11   Liu    A 1.00  1.42  2.18     12   Liu    B 0.67  0.76  1.64     13   Liu    C 1.39  1.62  3.38     14   Liu    D 0.56    NA    NA     15   Liu    E 1.00  1.58  2.02


title: Why would I want to use my data this way, and what am I supposed to do with this regression?

text: I'm writing an essay on the effects of M&A announcements on stock prices. I understand most of the requirements perfectly fine, but I'm stumped on one particular part of the essay. I have for every target and acquiror stock prices in the range [-170, 10] days, where day 0 is the announcement date. I use [-170, -71] to calculate the abnormal returns in [-10, 10] which is the event window I'm interested in. The next part confuses me: The idea underlying the in-between period [-70, -11] is that for each firm/stock 4 variables can be con-structed, that is the stock-price runup = stock return [-70, -11], the stock-sigma = stock vola-tility [-70, -11], the market-price runup = market return [-70, -11], and the market-sigma = market volatility [-70, -11]. Following this, I'm supposed to use these four variables for both univariate and multivariate regressions for CAR[K,L] where [K,L] will be chosen depending on the number of days for which a significant announcement effect has been found. My question is: why would I want to have this in my essay, I have literally no idea why this would be interesting. What is the purpose of the four variables, and why can the be used to regress on the CAR? I know this is quite an economical question, but I think the theory is statistical in nature, and there is not an economics stackexchange-site, so I decided to ask this question here. If anyone could provide any assistance, even the smallest hint, I would be most grateful. Thank you for any help you can give!


title: Why feature maps are indexed by two indices?

text: I'm reading about convolutional neural networks. As I understood a feature map is a set of neurons (i.e like a single hidden layer in traditional ANN). So why feature maps are indexed by (i,j)? Shouldn't a feature map be indexed only by one index because it's just a set of neurons? In the following picture layer C1 has 6 feature maps. The picture says that a feature map is 28 by 28. But why?! Why 2D?! ![enter image description here](http://i.stack.imgur.com/V79IH.png)


title: ordered probit in JAGS

text: I am running into an error               Error in node Y1[5]      Observed node inconsistent with unobserved parents at initialization      trying the following bugs model in JAGS accessed from R: Apparently `Y1[5] = 2` causes problems in `dcat`. But I can't seem to figure out why. Any help kindly appreciated. ### DATA               d1 <- matrix(c(1,1,1,1,2,2,2,3,3,3,3,3,4,4,5,5,5,                    0,0,0,0,0,0,1,0,0,1,1,1,1,1,0,1,1),ncol=2)     d2 <- matrix(c(1,1,1,2,2,3,3,3,4,4,5,5,5,                    0,0,0,0,0,0,0,1,1,1,1,1,1),ncol=2)     d3 <- matrix(c(1,1,1,1,1,2,2,2,2,3,3,3,                    0,0,0,0,0,0,0,1,1,1,1,1),ncol=2)     bugdata <- list(         N1 = 17,          N2 = 13,         N3 = 12,         ncat1 = 4,         ncat2 = 5,         ncat3 = 3,         Y1 = d1[,1],         Y2 = d2[,1],         Y3 = d3[,1],         D1 = d1[,2],         D2 = d2[,2],         D3 = d3[,2]     )      ### INITS               buginit <- list (              list (                 ta=1, tb=1,                  a=0.5, b=2, al2=0.1,                   alpha1 = 0.5, alpha2 = 0.5, alpha3 = 0.5,                  beta1 = 2, beta2 = 2, beta3 = 2,                  theta1_0=c(-.5,0,1),                  theta2_0=c(-.5,0,1,2),                  theta3_0=c(-.5,0)             )        )      ### MODEL               model{          # SET d1     for (i in 1:N1) {               for (j in 1:(ncat1-1)) {              probit(Q1[i,j]) <- (theta1[j] - beta1*D1[i])/exp(al2 + alpha1*D1[i])            }              p1[i,1]     <- Q1[i,1]          for (r in 2:(ncat1-1)) {              p1[i,r] <- Q1[i,r] -  Q1[i,(r-1)]         }          p1[i,ncat1] <- 1 - Q1[i,(ncat1-1)]               Y1[i] ~ dcat(p1[i,1:ncat1])          }           # SET d2     for (i in 1:N2) {               for (j in 1:(ncat2-1)) {              probit(Q2[i,j]) <- (theta2[j] - beta2*D2[i])/exp(al2 + alpha2*D2[i])            }         p2[i,ncat2] <- 1 - Q2[i,(ncat2-1)]          p2[i,1]     <- Q2[i,1]          for (r in 2:(ncat2-1)) {              p2[i,r] <- Q2[i,r] -  Q2[i,(r-1)]         }          Y2[i] ~ dcat(p2[i,1:ncat2])          }           # SET d3     for (i in 1:N3) {               for (j in 1:(ncat3-1)) {              probit(Q3[i,j]) <- (theta3[j] - beta3*D3[i])/exp(al2 + alpha3*D3[i])            }         p3[i,ncat3] <- 1 - Q3[i,(ncat3-1)]          p3[i,1]     <- Q3[i,1]          for (r in 2:(ncat3-1)) {              p3[i,r] <- Q3[i,r] -  Q3[i,(r-1)]         }          Y3[i] ~ dcat(p3[i,1:ncat3])          }           # PRIORS...          for (t in 1:3) {theta1_0[t] ~ dnorm(0, 0.1)}     theta1[1:3] <- sort(theta1_0)     for (t in 1:4) {theta2_0[t] ~ dnorm(0, 0.1)}     theta2[1:4] <- sort(theta2_0)     for (t in 1:2) {theta3_0[t] ~ dnorm(0, 0.1)}     theta3[1:2] <- sort(theta3_0)          alpha1 ~  dnorm(a,ta) T(-10,10)      alpha2 ~  dnorm(a,ta) T(-10,10)      alpha3 ~  dnorm(a,ta) T(-10,10)           beta1 ~  dnorm(b,tb)       beta2 ~  dnorm(b,tb)       beta3 ~  dnorm(b,tb)            al2  ~ dnorm(0.0, 0.1) T(-10,10)      a    ~ dnorm(0.0, 0.1) T(-10,10)      b    ~ dnorm(0.0, 0.1)      ta   ~ dgamma(0.01, 0.01)      tb   ~ dgamma(0.01, 0.01)           }


title: Best subsets with logistic regression in R?

text: @Dmitrij Celov posted an answer to the question How to do logistic regression subset selection saying that this paper was a good reference. I was wondering if anyone knows if any of these authors ever created an R package based on it.


title: Two-way ANOVA robustness against normality violations

text: I recently read the article "Is It Really Robust?" (Schmider, E., Ziegler M., Danay, E., Beyer, L. & Bühner M., 2010) The authors of the article come to the conclusion that an ANOVA is quite robust against violations of the normality assumption. However the study focusses on one-way ANOVAs. I would like to know if these findings can be generalized on a two-way repeated measures ANOVA? References to peer-reviewed articles are highly appreciated! * * * ### References Schmider, E., Ziegler M., Danay, E., Beyer, L. & Bühner M. (2010). Is It Really Robust?. _Methodology, 6 (4),_ 147 - 151.


title: Statistical test for finding significant positions having deviated values

text: I have near about 50 files(each file corresponds to a patient) with 4 columns -               chromosome     start.position       stop.pos       value First 3 columns in all 50 files are same and fourth column is value which is different (or can be same for some) for all patients. Basically this value corresponds to copy number (structural variation in genome). I want to apply some test (or if you can suggest any other procedure) to find out values which are deviating as compared to others. Example-               file1               chromosome     start.position       stop.pos           value           1              10                  110             4           2              100                 200             5           2              500                 600             0.5       file2               chromosome     start.position       stop.pos           value           1              10                  110             2.5           2              100                 200             6           2              500                 600             0.6       file3               chromosome     start.position       stop.pos           value           1              10                  110             3           2              100                 200             5.5           2              500                 600             3.5       file4               chromosome     start.position       stop.pos           value           1              10                  110             2           2              100                 200             0.9           2              500                 600             3       file5               chromosome     start.position       stop.pos           value           1              10                  110             8           2              100                 200             4.5           2              500                 600             2.5       So in the output, file5 has significant value (as compared to same position of other files) at row1, file4 has significant value at row2 and file1 and file2 has significant value at row3. I have 1 solution in my mind- pick a row from each file and take the average of values. Then compare this average with each file and return the position if it is significant as compare to average. But I feel there should be some better solution for this. Can you please suggest any test or any other procedure that you will use to complete this task. Thanks in advance. EDIT: After reading Joel's comment. Yes I have a specific question in my mind. As I mentioned these files corresponds to each patient, so I am trying to find if a patient has significantly increase(or decrease) in value at certain position (i:e row). Ok we can talk about only 1 row now. Just consider 1 position (1 row), take the values ( in this case, for first row we will have 5 values as we have five files) and from these values check which value is significantly deviating from others. So if you can suggest me solution for 1 row, I can do it for all the rows one by one. Please let me know if it is not clear.


title: How to find P(X > F) in ANOVA F-test?

text: _(Apologies for the ASCII tables, Stackexchange doesn't allow HTML tables and since I'm not supposed to link to an image, this is the only way I know of showing the data)_ I'm learning about ANOVA F-testing and stumbled upon this problem: There exists the following set of data regarding four teaching methods and the scores of students who were subject to each teaching method:               Method 1  Method 2  Method 3  Method 4     ----------------------------------------       65        75        59        94       87        69        78        89       73        83        78        80       79        81        62        88       81        72        83       69        79        76                 90     ------------------------------------      $\bar{x} = 75.67$       $78.43$          $70.83$         $87.75$   $s = 8.17$          $7.11$            $9.58$           $5.80$   Where $\bar{x}$ is the mean of all scores for each teaching method and $s$ is the standard deviation for each method. After conducting an ANOVA, I derive the following ANOVA table:               Source     | Deg. Freedom |    SS     |    MS     |   F  |     -------------------------------------------------------------     Treatment  |     3        | 712.59    |  237.53   | 3.77 |     Error      |     19       | 1196.63   |  62.98    |  -   |     Total      |     22       | 1909.22   |    -      |  -   |      Now the question is: Test a level $\alpha = 0.05$ The null hypothesis is that there is no difference in mean achievement for the four teaching techniques. So to restate: $H_0$: Teaching technique _does not_ have an influence on mean achievement of students   $H_1$: Teaching technique _does_ have an influence I work out the critical f value to be $f_{3,19;0.95}$ from an F-distribution table to be $3.1274$ The next step is what I don't understand. We can claim that the teaching technique does have an influence on the mean achievement of the students (with less than 5% chance of being wrong)   The associated p-value is: $p = P(X>3.77) = 0.0281$ and indeed, p < $0.05$ (hence reject of $H_0$) But now where is this $0.0281$ from? It looks to be the probability that X > 3.77. If I'm not wrong, in this case $X ~ F_{3,19} = 3.1274$ (as calculated before) and so the p value should be the probability that 3.1274 is greater than 3.77. Now how can 3.1274 ever be greater than 3.77?


title: Uniform Random on $(-\infty,\infty)$

text: Imagine picking a 1 when any real number is equally likely. What is the pdf? Does this idea have a known use? What is its name? There could be a use for a uniform random real number. It could end up always being integrated. It could be a prior for a number that is completely unknown. It would be like a delta function, which integrates to a finite number though it doesn't have a well-defined value at zero. The delta function is greater than any number at zero. This uniform random on (-inf,inf) would be not well- defined anywhere. It would be less than any positive number and greater than zero. What is it?


title: What  causes loss of power in hypothesis testing?

text: I am just wondering what kind of things could cause power to be lost in general (the opposite of reducing the noise of the data to increase the power of the test). Could anyone give some summary or ideas how to think of the power issue?


title: How to statistically analyze significance with 7-point Likert scale data?

text: I am working on a business research paper. The hypothesis is that the sales rep influences the customer's decisions. Independent variables consist of the sales rep's: knowledge, pushiness, friendliness, etc. Dependent variable is the customer's willingess to buy. In all, I have 12 questions that I asked on a Likert scale ranging from 1 (not important at all) to 7 (extremely important). At this point, I have all the means and standard deviations calculated but am not sure how to proceed to evaluate whether the whole set of independent variables proves my hypothesis that the rep's qualities does influence the customer to buy. I've done some reading on regression, ANOVA and chi squared but would really like to understand how I can apply them to this situation. Additionally, I want to know if I compute the composite score of the subgroups, i.e. men versus women, how can I compare them statistically? Thanks for any and all suggestions!


title: Model function for discovering irrelevant dimensions with L1 regularization

text: For homework I have been given a 20-dimensional input $x \in \mathbb{R}^{20}$, many of which are suspected to be irrelevant. I tried using L1-norm Lasso regularization to uncover which dimensions contribute to the output: $$L(\beta) = \sum_{i=1}^n (y_i - \phi(x_i)^T \cdot \beta)^2 + \lambda \sum_{j = 1}^k l(\beta_j)$$ Please note, that instead of $|\beta_j|$ another function is used, where $$l(\beta_j) = \begin{cases} |\beta_j| - \varepsilon/2 & \textbf{if } |\beta_j| \geq \varepsilon\\ |\beta_j^2| / (2\varepsilon) & \textbf{if } |\beta_j| < \varepsilon\\ \end{cases}$$ With the resulting gradient: $$\frac{\partial}{\partial\beta} L(\beta) = -2 \sum_{i=1}^n - \phi(x_i) \cdot (y_i - \phi(x_i)^T \cdot \beta) + \frac{\partial}{\partial\beta_m}\sum_{j=1}^kl(\beta_j)$$ $$ \frac{\partial}{\partial\beta_m}\sum_{j=1}^kl(\beta_j)= \begin{cases} 0&m\gt k\;,\\ sign(\beta_m)&m\le k\;,|\beta_m|\ge\epsilon\;,\\ sign(\beta_m)/\epsilon&m\le k\;,|\beta_m|\lt\epsilon\;. \end{cases} $$ In order to find the minimum I applied gradient descent on the differentiated Lasso function and received a $\beta$ vector after 6 to 10 iterations. However, I don't see how this helps me to uncover the irrelevant dimensions. How should I proceed?


title: Estimating parameters in a model with a periodic design parameter

text: We recently studied a model with a likelihood function of the form $$\Pr(d_i|\omega,t_i)=(1-d_i)+(2d_i-1)\cos^2(\omega t_i),$$ where $d_i\in\{0,1\}$, $\omega$ is an estimation parameter, and $t_i$ are determined by the experimenter. This is a set of independent but non- identically distributed Bernoulli trials (think: trying to estimate the time- varying bias of a coin). Since this was a simple one estimation parameter, one design parameter model, we could essentially employ a first principles approach to the theory and computations (the paper is here: http://arxiv.org/abs/1110.3067). Now that we plan to generalize this to multiple estimation and design parameters, I really want to make sure we are doing the most efficient thing possible. I did a bit of searching but couldn't find anything in the literature on statistical models that depend on a periodically varying design parameter $t$. Has anyone seen a model like this studied before?


title: what are the dependent and independent factors and categorical for below project

text: my project is identification of factors influencing motor bikes for home to office trip. there are some attiributes like gender,age,marital status,education qualification,job type,working sector, no. of members in household,working members in household,monthly income,total household income,no. of motor bikes,members having motor bike driving licence(D.L),no. of cars,members having car D.L,do you have D.L for motor bike, do you have D.L for car,type of house, live-in status of house, do you work from home, do you use internet for office purpose, no. of days you work from house,purpose of trip, timings of office i.e., starting and ending time, vehicle for above trip,distance from home to office,time from home to office and office to home,reason for preferring to private vehicle,likely fare (if you go by bus or shared auto),drop off any kid to school while going to office,does the school come on your way to office, spouse is working person,do you drop off your spouse to her office,city bus is accessible to your house,walking distance from home to bus stop,waiting time at bus stop... and some of these questions are related to one another like spouse is working ?? if no he should not fillup the next question. also like drop off any kid to school.. wat should i use in the related question...wheather '-' (or) 9999...also i am totally confused about the dependent and independent and categorical variables........ also wat the model we are going to use for this project wheather binary logistic (or) multinomial logistic ...


title: Derivation of closed form lasso solution

text: For the lasso problem $\min_\beta (Y-X\beta)^T(Y-X\beta)$ such that $\|\beta\|_1 \leq t$. I often see the soft-thresholding result $$ \beta_j^{\text{lasso}}= \mathrm{sgn}(\beta^{\text{LS}}_j)(|\beta_j^{\text{LS}}|-\gamma)^+ $$ for the orthonormal $X$ case. It is claimed that the solution can be "easily shown" to be such, but I've never seen a worked solution. Has anyone seen one or perhaps has done the derivation?