title: Consequences of an improper link function in N alternative forced choice procedures (e.g. 2AFC)?

text: Background: In some cognitive psychology research areas N-alternative forced choice tasks are common. The most common of these is a two alternative forced choice (2AFC). This usually takes the form of participants being given a stimulus and asked to make one of two judgement, e.g. the target stimuli is present/absent, the stimulus on the left is the same/different than the the one on the right, etc. Designs in which the experimental data is from a 2AFC but there is only one data point per subject are rare, but do exist, e.g. some eye-witness identification research. Since the dependent variable (correct/incorrect) is binary, these experiments are reasonable places to use logistic regression. My question is this: since chance performance is 50% in a 2AFC trial, is it still reasonable to use the standard logistic link function? Specifically, the logistic function has a minimum value approaching 0% correct, but in practice participants in a 2AFC should be correct at least 50% of the time due to chance. I imagine the following case in which it may present a problem: an independent variable is assessing the difficulty of the discrimination (e.g. difficulty 1, easy - 5, hard; please note this is introduced in ordinal terms only for ease of comprehension - for the sake of this discussion consider this variable as being interval) - participants got 50% correct at 5 and 4, 75% correct at 3, 85% correct at 2, and 99% correct at 1. Would using a standard logistic link function cause us to underestimate the slope? [I think so, but please correct me if I'm wrong, see below] Edit: Those who have answered my question so far have expressed that the way in which I set up the problem was unclear. I'm providing the sample below to help clear things up.               library(psyphy)     make.data <- function(zero,one)         {             return(c(rep(0,zero),rep(1,one)))         }     center <- function(x) {return(scale(x,scale=FALSE))}     logit.data <- data.frame(Score=c(make.data(50,50),make.data(50,50),make.data(25,75),make.data(15,85),make.data(1,99)), Difficulty=rep(5:1,each=100))     logit.data$Difficulty2 <- center(logit.data$Difficulty)^2     standard <- glm(Score~center(Difficulty),data=logit.data,family=binomial) #standard link function     standard.2 <- glm(Score~center(Difficulty)+Difficulty2,data=logit.data,family=binomial) #standard link function, but better with a quadradic     revised.link <- glm(Score~center(Difficulty),data=logit.data,family=binomial(mafc.logit(2)))     AIC(base)     AIC(base.2)     AIC(revised.link)     coef(base)     coef(base.2)     coef(revised.link)     #plot     plot(diffs,plogis(coef(standard)[1] +coef(standard)[2]*center(diffs)),xlab="Difficulty",ylab="Pr(Correct)",ylim=c(0,1),col="blue",type="l");abline(.5,0,col="Orange");lines(diffs,plogis(coef(standard.2)[1]+coef(standard.2)[2]*center(diffs)+coef(standard.2)[3]*center(diffs)^2),col="Cyan");lines(diffs,(p2afc(coef(revised.link)[1]+coef(revised.link)[2]*center(diffs))),col="Green");lines(5:1,c(.55,.60,.75,.85,.99),col="Black")      ![alt text](http://i.stack.imgur.com/elxVG.jpg) In the above image the orange horizontal line marks 50% correct responses. The jagged black line represents the data supplied to the estimation equation (note the values for 4 and 5 disappear behind the orange 50% marker). The blue line is the equation produced by a standard logistic link. Note that it estimates below 50% accuracy when discrimination is most difficult (5). The cyan line is the standard logistic link with a quadratic term. The green line is a non-standard link that takes into account that the data comes from a 2AFC experiment where performance is very unlikely to fall below 50%. Note that the AIC for a model fit using a non-standard link function is superior to the standard logistic link function. Also note that the slope for the standard equation is less than the slope for the standard equation with the quadratic term (which more accurately reflects the real data). Thus, using a logistic function blindly on 2AFC data does (at least) appear to underestimate the slope. Is there a problem with my demonstration that means that I am not seeing what I think I am seeing? If I'm correct, then what other consequences (if any) are there of using the generic logistic function with 2AFC data [presumably extensible to NAFC cases]?


title: Mixed regression models and custom link functions in R?

text: It seems like the current revision of lmer does not allow for custom link functions.   1. If one needs to fit a logistic linear mixed effect model with a custom link function what options are available in R?   2. If none - what options are available in other statistics/programming packages?   3. Are there conceptual reasons lmer does not have custom link functions, or are the constraints purely pragmatic/programmatic?


title: Computing the Gini index

text: ![enter image description here](http://i.stack.imgur.com/CLm9y.jpg) How do I compute the Gini index using **Instance** attribute as attribute test condition? I calculated the Gini, but I have no clue how to do it for this Instance attribute. $$\text{Gini for } a_1 = 0.345 $$ $$\text{Gini for } a_2 = 0.493 $$ $$\text{Gini for } a_3 = ?$$ I am guessing the answer to this is that Instance attribute has no information gain. However, I can't prove this. $$\text{Gini} = 1 - \sum_i p(i|t)^2$$


title: R - show movement of points

text: Suppose I have data of the following form               t1 = data.frame(a=c(2,3,3,1,5),b=c(6,4,5,2,1))     t2 = data.frame(a=c(3,4,4,1,8),b=c(5,5,5,3,3))      If you consider `plot(t1$a,t1$b)` and `plot(t2$a,t2$b)`, you can imagine that the second plot is produced by taking every point on the first plot and moving it to its new location in the second plot. I'd like some nice way of visualizing this. One way I thought of was to represent each `t1` `(a,b)` as a vector reaching to its corresponding `t2` `(a,b)` point, and then plotting a bunch of vectors, sort of like what I used to do in differential equations (I'm forgetting the technical name for this type of plot). Basically, it would be a plot of a bunch of little arrows that originate from each `t1` `(a,b)` and point to each `t2` `(a,b)`, respectively. Any ideas of how to do this in `R`, or any other ideas for visualizing and looking at how each dot moves. NOTE: in the real version of what I hope to accomplish, I have lots of data, maybe 200 points, and I'm trying to show that in general, points are moving towards smaller `a` values and smaller `b` values from `t1` to `t2`. But I really want to show this some other way besides something like a histogram of the difference between `t2` and `t1` average of `a` and `b`.


title: Estimating population size of a subgroup based on independent samples without replacement

text: Let a bag have 1000 balls, 100 red and 900 blue. Now, let us get ten independent samples (without replacement within each sample) of the original population (10% each). Proceed to count the number of red balls in each sample, i.e.:               (9, 10, 10, 11, 13, 8, 5, 15, 12, 9)      To estimate the number of red balls in the original population, one can calculate the average of the above list and divide by the sample. In this case:               (9 + 10 + 10 + 11 + 13 + 8 + 5 + 15 + 12 + 9) / 10 / 0.1 = 102      My questions are thus:   1. Instead of a number, how can I obtain the probability mass function of the estimated population? Is it sound to use the central limit theorem here, where the distribution $N = (\mu,\sigma^2/n)$ ? I have a problem with this approach since negative values don't make sense.    2. Now imagine I have an infinite numbers of samples (all w/ 10% of the population), or at least a very high (but known) number of samples, and keep getting the result for each of them, one by one. How do I decide when it is likely ok to stop (within a certain confidence)? * * * Some clarifications over the original problem:   1. The population size is known _à priori_.   2. But the number of colors is not known _à priori_.    3. Neither their ratio.   4. Each sample take balls without replacement.   5. Samples are independent (after a sample is taken, all balls are replaced).


title: How to determine whether a non-stationary time series has change-in-mean or change-in-variance?

text: Conditions of a weak-stationary time series are: (1) constant in mean (2) constant in second-order moments. We can test whether a time series is stationary using like `kpss.test()` in R. If say, a time series $X_t$ is non-stationary by the test. Is there a way to check whether it is not constant in mean, or, its second-order moments are changing over time?


title: Cumulative Distribution Function Inequality (Discrete Distributions)

text: Let a discrete Random Variable $T$ have CDF $F_T(T)$. Could you please help me understand why $$ P \left[ F_T (T) \leq a_1 \right] \leq a_1 $$ I know that the result holds with equality for the continuous case, it is known as Probability Integral Transform, but I am having trouble understanding it for the Discrete one since the inverse is not defined. Thank you.


title: R libraries for deep learning

text: I was wondering if there's any good R libraries out there for deep learning neural networks? I know there's the `nnet`, `neuralnet`, and `RSNNS`, but none of these seem to implement deep learning methods. I'm especially interested in unsupervised followed by supervised learning, and using dropout to prevent co-adaptation.


title: Does it matter if I use correlations or regression coefficients to suggest areas to focus on to improve overall customer satisfaction?

text: I am working with customer satisfaction data where the dependent variable is "Overall satisfaction" and the independent variables are satisfaction with various areas such as customer support, delivery etc. I want to suggest areas where the company should focus on in order to improve overall satisfaction. Option 1: I could look at correlations between the 'Overall satisfaction" and the independent variables and suggest that the company focus on the top 3 positive correlations as areas for improvement. Option 2: I can use a linear regression and suggest that the company should focus on the areas associated with the 3 highest regression coefficients. Are the two options equivalent? If not, which one is the better approach?


title: Similarity of new element x with the training set X

text: Suppose we have trained a model (function, algorithm) $M$ which gives prediction to a new sample $x$ not observed in the training set, $M(x)$. It is natural to assume that the quality of prediction $M(x)$ depends on similarity of $x$ with the training set $X$. I wonder if there are any indices/methods which would evaluate similarity of $x$ with training set $X$, or any indices which would evaluate the quality of prediction $M(x)$. To be more specific, I don't want an overall performance of the model (which could be evaluated by RMSE, R-squared or deviance function), but I would like to estimate the quality of the prediction $M(x)$ for this particular $x$. It would be great if anyone could point me to the literature on this problem.


title: Sample size and correct choice of test in g*power

text: I am new to g*power and have a question about which test I have to choose and how to interpret the given sample size. I have 2 measurements (pre / post), one control-group and one intervention- group. In a reference study, I was able to find an effect size of .7 for the most important parameter. In g*power I chose F-test -> ANOVA RM within-between interaction, a power of 0.8 and alpha .05. Here are my questions:   1. Did I choose the correct test?   2. Is the resulting sample size referring to each group or total?   3. Regarding the effect size of .7 I found in a reference study, am I allowed to fill it in g*power -> effect size f, or should I leave it as it was at .25?


title: How can I calculate the statistical power of a MI based likelihood ratio test with many df?

text: I'd like to calculate the statistical power of a likelihood ratio test based on mutual information. I'm using the asymptotic property that the test statistic $G = 2NI$ where $N$ is the sample size and $I$ is the mutual information. I have a large sample size (say $10^6$ observations) but also many degrees of freedom (say $6000^2$), so I am calculating the critical value of the $\chi^2$ dist using the normal approximation $\chi^2 \sim \mathcal{N}(\nu, 2\nu)$. I suspect I am trying to show that the required sample size for a significant confidence with so many degree of freedom is unreasonable so any alternative approaches to argue this point also welcome. I'm not a statistician so apologies if this question is poorly phrased or misthought!


title: What are good econometrics handbooks?

text: There's a question for textbooks here, but I would like to ask similar question about handbooks: What econometrics handbooks would you recommend? Assumed audience is researchers and graduate-level students. It needs to include the material of Greene's _Econometric Analysis_ and Wooldridge's _Econometric Analysis of Cross Section and Panel Data_ in a denser form (without discussions and proofs), as well as mainstream techniques not mentioned in these two books. One option is Elsevier's 6 volumes of _Handbook of Econometrics_ series edited by Griliches and Instriligator. However, perhaps you would recommend other handbooks, more concise or otherwise? I suggest two possible formats: one is a reference-card format with minimum explanations, and the other one is a more extended format with proofs and more detailed exposition.


title: determining sample size from a pooled sample

text: Let's say you have an overall group of financial institutions (fi) made of banks (ba) and federal credit unions (fcu). The total number of fi is 73. The number of ba is 46 and the number fcu is 27. Let's assume also the total numbers of employees working for ba is 2500 and the total numbers of employees working for fcu is 1080. How will you construct a sample or sample groups that are representative of each group; that is, ba and fcu?


title: How does one determine what ARL0 should be used on CPM package to test for Structural Change

text: I'm trying to find multiple break points by using `processStream` from `CPM` package on R. Can someone enlighten me on what is `ARL0` how does one determine what `ARL0` should be used for?               processStream(ret.fin.chn,"Kolmogorov-Smirnov",ARL0=500,lambda=NA)     $changePoints         [1]   59   75  250  286  443  448  663 1037 1042 1261 1576 1842 1853 2013 2035 2621 2633         $detectionTimes     [1]   73   89  285  334  447  503  670 1040 1145 1428 1639 1951 1874 2030 2078 2632 2644      while               processStream(ret.fin.chn,"Kolmogorov-Smirnov",ARL0=2000,lambda=NA)     $changePoints         [1]   59   75  663 1037 1261 1559 1842 2013 2035 2621 2633         $detectionTimes     [1]   75   90  691 1041 1480 1688 2026 2032 2266 2633 2646      and               processStream(ret.fin.chn,"Kolmogorov-Smirnov",ARL0=3000,lambda=NA)     $changePoints         [1]   59   75  663 1037 1261 1559 1842 2013 2149         $detectionTimes     [1]   75   92  692 1041 1490 1690 2026 2032 2284      Tt seems that different `ARL0` will give fewer break point detection, is that a good thing? **Note** : the time series `ret.fin.chn` contains 2749 rows. Below are excerpts from `R` help: > ARL0 > > Determines the `ARL_0` which the CPM should have, which corresponds to the > average number of observations before a false positive occurs, assuming that > the sequence does not undergo a chang. Because the thresholds of the CPM are > computationally expensive to estimate, the package contains pre-computed > values of the thresholds corresponding to several common values of the > `ARL_0`. This means that only certain values for the `ARL_0` are allowed. > Specifically, the `ARL_0` must have one of the following values: 370, 500, > 600, 700, ..., 1000, 2000, 3000, ..., 10000, 20000, ..., 50000


title: Weighted Least Squares Estimate

text: Here is a problem from a practice test. Suppose that $$X_i = \mu + \epsilon_i,\quad i=1,\ldots,n\quad \epsilon_i\sim N(0,\sigma^2_1)$$ $$Y_i = \mu + \delta_i,\quad i=1,\ldots,m\quad \delta_i\sim N(0,\sigma^2_2)$$ All $\epsilon_i$'s and $\delta_i$'s are independent. The paramters $\mu, \sigma_1^2, $ and $\sigma_2^2$ are unknown. Let $\theta=m/n$, $\rho=\sigma_2^2/\sigma_1^2$. Suppose $\rho$ is known. Show that the least squares (weighted) estimator of $\mu$ is $$ \hat{\mu} = \dfrac{\rho\bar{X} + \theta\bar{Y}}{\rho+\theta}$$ MY ATTEMPT: I can't figure out how to use the fact that $\rho$ is known. I tried $$\hat{\mu} = \text{argmin}\left\{\sum_{i=1}^n (X_i-\mu)^2 + \sum_{i=1}^m (Y_i-\mu)^2\right\}$$ and arrived that the weighted averaged $$ \hat{\mu} = \dfrac{n\bar{X} + m\bar{Y}}{n+m}$$ But again this does not use the fact that we know what the ratio $\sigma_2^2/\sigma_1^2$ is. Any ideas?


title: How to express joint conditional probability with multiple conditions

text: I want to express the joint probability of $\Phi_A$ and $\Phi_B$: $p(\Phi_A, \Phi_B)$ conditioned that $\Phi_A$ and $\Phi_B$ are both greater than some value C. How would I express this mathematically? I guess my intuition says: $p(\Phi_A, \Phi_B | \bf{\Phi} >C)$ Is this correct? Is there a better way to express this?


title: Failure analysis on server

text: If you could review and share any insights/guidance you may have to share, i would sincerely appreciate it!!!!! **Dataset** Server id---------Repair id---------Repair Type---------Date of repair   * each server can have the same type of repair multiple times in its lifetime    * server id identifies a server    * repair type determines model of part being repaired ( disk model/ram model )   * there are around 600000 servers with repair data exceeding 3 million records over 3 years **Goal**   * Determine repair types that may be causing more frequent repairs subsequently  **Constraints** I feel there is no certain way of establishing causality - hence let the business know - but i do feel there should be a way of determining the failure types that are followed by a higher frequency of failures **My current approach** 1) For every failure type for every server, determine                     R = X / Y                Where                  X is              ( number of failures after the failure in question until the last failure  )              ____________________________________________________________________________              ( number of days after the failure in question to the last failure )                  Default to 0 if there are no failures                       Y is               ( total number of failures until the last failure  )              ____________________________________________________              ( number of days from server start to the last failure )                 Default to 0 for no failures       2) Retain only the first value of R for each failure type per server 3) Determine the average and variance for the R values of each failure type 4) Sort by average desc, variance asc and filter the first N failure types **Your guidance: ???**


title: Pinsker's Inequality for Bayesian hypothesis testing

text: I am wondering if one can relate the KL divergence to the probability of error in a Bayesian binary hypothesis testing setting. That is, we have to decide between hypotheses $A$ and $B$ given observationx, and we have prior probabilities of either being true $p(A)$ and $p(B)$ such that $p(A)+p(B)=1$. Under either $A$ or $B$, observation induces (different) probability distributions $p_A(x)$ and $p_B(x)$. The setting is identical to a previous question, and this question may be considered its extension. The following restatement of eq. (3) from the previous question lower-bounds the probability of error of the hypothesis test: $$P_e = \frac{1-||p_A(x)p(A) - p_B(x)p(B)||_1}2\tag{1}$$ where $||\cdot||_1$ denotes the $L_1$ norm. When the prior probabilities are equal, that is $P(A)=P(B)=\frac{1}{2}$, (1) reduces to the expression involving the total variation distance between probability distributions induces by the two hypothesis: $$P_e = \frac{1-\frac{1}{2}||p_A(x) - p_B(x)||_1}2\tag{2}$$ Pinsker's inequality states that $$||p_A(x) - p_B(x)||_1\leq\sqrt{2D(p_A(x)||p_B(x))}\tag{3}$$ where $D(p_A(x)||p_B(x))=\sum_x p_A(x)\log\frac{p_A(x)}{p_B(x)}$. Of course (3) can be used to lower-bound (2) as follows: $$P_e \geq \frac{1-\frac{1}{\sqrt{2}}\sqrt{D(p_A(x)||p_B(x))}}2\tag{4}$$ However, I am wondering if a similar KL-divergence based lower bound can be derived for (1). Does weighting by the priors work in this case, i.e. is (5) true: $$||p_A(x)p(A) - p_B(x)p(B)||_1\leq\sqrt{2D(p_A(x)p(A)||p_B(x)p(B))}\tag{5}$$ I don't think it is true. So, is there a KL-divergence based lower-bound for (1)?


title: Data set of means

text: Suppose we have a data set of means and standard deviations. So the data are of the form $X_1=(\bar{x}_1, \bar{x}_2,\dots, \bar{x}_n)$ and $X_2 =(s_1,s_2, \dots, s_n)$. Is there any utility in looking at $\overline{X_1}$? Also is there a relationship between $\text{Stdev}(X_1)$ and $X_2$? It seems that $X_2$ gives intervariability and $\text{Stdev}(X_1)$ gives the intra- variability. Maybe this is related to the total law of variance?


title: How to determine post-treatment weight with uneven time points?

text: I am looking for suggestions on how to deal with uneven time points and missing data in a single-group repeated-measures design. **BACKGROUND:** A group of 220 subjects went through a 4-month weight loss treatment, and were given wireless scales to weigh themselves daily at home. The primary outcome is change in weight from pre- to post-treatment (for which I would normally use a paired samples t-test on). **HOW TO DETERMINE POST-TREATMENT WEIGHT?**   We have pre-treatment weights for everyone, but post-treatment weight is hard to determine since not all the subjects weighed-in at the same time after treatment (and some dropped out altogether). **APPROACH A: Creating "Time Windows" to Determine Post-Treatment Weight.**   I've tried creating arbitrary time windows (1, 2, 4, or 6 weeks around the treatment end date), and then taking the median weight to determine "post- treatment weight." I found (all p's < 0.01):   1 week window: $n = 124$, $-9.55$% weight loss   2-week window: $n = 143$, $-7.52$% weight loss   4-week window: $n = 161$, $-6.39$% weight loss   6-week window: $n = 168$, $-6.30$% weight loss   However, this is problematic because it excludes dropouts. Furthermore, the bigger the window, the bigger the "n", but smaller the effect size (because subjects who did not weigh in regularly also lost the least weight) and less valid the weight in terms of being measured at the same time between subjects. **APPROACH B: Intention-to-Treat with Last Observation Carried Forward.** Just using the last observation of all the subjects, up to 2 weeks post-treatment:   $n = 220$, $-4.22$% weight loss However, this is problematic because the dropouts are depressing the effect size (though retains the largest "n"). **APPROACH C: MIXED MODELS?**   I've read that mixed models can handle uneven time points and missing data, though I'm not super familiar with this procedure, and am not sure what are the fixed and random are effects here. My attempts to run it in SPSS 21 caused an "insufficient memory" error. **APPROACH D: SOMETHING ELSE? (Multiple Imputation/Regression/ANCOVA, etc.)**


title: Covariance of a bivariate Gaussian given identity matrix

text: I have a homework problem about finding an optimal decision boundary. I know the formula (not really the process) for calculating one, so that may be another question entirely, but I do know I need the mean, covariance, and priors. The question is shown below               Suppose points in R^2 are being obtained from two classes, C1 and C2, both of which       are well described by bivariate Gaussians with means at (0,0) and (1,3) and covariances      I and 2I respectively. I is the (2x2) identity matrix. If the priors of C1 and C2 are      0.4 and 0.6 respectively, what is the ideal (i.e. Bayes Optimal) decision boundary      (derive the equation for this boundary)?      I'm puzzled when it says "covariances I and 2I... I is the (2x2) identity matrix." Is there a way to calculate a numerical value for covariance based on this information? If so, how? I've never heard it phrased this way.


title: Sensitivity analysis for modeled time series

text: I have very basic knowledge of stats so my question may sound very simplistic. I have a large time series of measured data and have calibrated a model with five parameters to make predictions. I would like to test the sensitivity of my model but I am not sure what statistical method to use. I've been reading about sensitivity analysis, however the output of my model is a time series with fluctuations and a simple sensitivity analysis does not seem appropriate since it looks at individual output points and not the entire series. Any recommendations on methods or approaches?


title: How do you treat 'don't know' answers for a composite likert scale?

text: I'd like to create a composite item from variables with the following scale: strongly disagree, disagree, neither agree or disagree, agree, strongly agree and don't know. I was going to code the scale as 1=strongly disagree, 5=strongly agree, and calculate the composite item with these values. However I don't know how I would code 'don't know' answers (which are about 2% of total answers in the overall sample for all variables). Any suggestions on how to do this or should I just be treating them as missing data and not include them in my composite item?


title: How to correct for correlation at baseline between predictor and "DV"?

text: So say I am testing the predictive value of a predictor towards a "DV" within a survey data-set, between two time-points. Say also that I have measurements of both variables at all time points and that I can see that the DV and the predictor correlate at baseline. If I use linear regression or pearson correlation to evaluate the associaton between the predictor at baseline and the DV at a later time-point, do I need to correct for the correlation at baseline? I figure that if I do not correct, somehow, for this initial correlation I am not only testing prediction of the DV but also the test- retest-reliability of the predictor. Hope I am making sense.


title: What is better for time series prediction: AR or ARIMA?

text: I am trying to make a prediction in a time series with window 512 and horizon 2. I want to know if it's worth using ARIMA, that seems to be hard to understand, instead of the simple Autoregressive model?


title: How to interpret R output for two-way repeated measures ANOVA

text: How do I interpret the following output from R for a two-way repeated measures ANOVA               Model <- aov(value ~ BWRatio * NumGraBoolean + Error(participant/(BWRatio * NumGraBoolean)))      ![enter image description here](http://i.stack.imgur.com/kQCo7.png) In particular, what does all these Error terms mean and what does their p-value signify? Please excuse my lack of understanding repeated measures and ANOVA in general, but what do the residuals here represent as well?


title: Determining order of AR model

text: Suppose that we have the following model $$ y[t] = A_1\sin(\omega_1 t+\phi_1)+A_2\sin(\omega_2 t + \phi_2)+ \cdots + A_p \sin(\omega_pt + \phi_p) + z(t) . $$ Let us call this signal as B. Then in matlab I executed the following commands               >> [pxxr,fr]=pyulear(B,50,1024,100);     >> plot(fr,pxxr)      I get the following picture ![enter image description here](http://i.stack.imgur.com/omDoL.png) It seems that order is correct, but I have done it using partial autocorrelation ploting. I know that there are other criterias for determining AR order, like AIC, BIC, Minimum discription length (MDL), but how to calculate them in matlab? Please help me with code or by some other criterias in matlab to determine order of AR model.


title: Which steps have to be done before fitting logistic curve to time-series?

text: I want to cluster time-series concerning sales of products. In the database I have 26weeks after launching each products and units sold each week. One of the method of clustering is to cluster parameters of growth curve fitted to time-series. Can somebody tell me what should I do before fitting logistic curve to time- series? Should I use standardization, create model to each time series, use lsm to count parameters and then cluster them? I would be grateful for any explanation.


title: Posterior Density in R

text: I'm new to the site, and to Bayesian statistics and was hoping to get some help. I'm currently working through some study exercises and am required to compute the mean and variance of the posterior density in R. The prior density is represented by a normal curve with a mean of 6 and variance of .25, whereas observed mean from a sample (from a normal density) of n = 24 is 7.69 with variance of 2.0. Most of the readings and study I've done relate to probability estimates - nothing specifically related to the exercises I'm currently trying to work through. Any assistance or advice or information would be greatly appreciated. Thanks


title: How to get the regression from a plot?

text: I have a dependent variable `C` and an independent variable `VPT`. `VPT` is the average volume per tree (in cubic foot) of a timber stand. `C` ist costs resulting from a cost function, that calculates the harvest costs per cubic foot (in $). The plot between the two looks like this: ![enter image description here](http://i.stack.imgur.com/LqfIs.png) The abrupt change at roughly 28 cubic feet results from interaction within the cost function. I want to create a regression model for the two. Based on the graph I assume it must be something like this:               fit <- lm(C ~ I(VPT^(-1)))      Is there a way to determine what the best fit is? Do I just have to try around until I get something that fits, or is there a more efficient way? * * * The cost function takes in four parameters:   * Volume per tree (VPT)    * Trees per acre (TPA)    * Skidding Distance (SD)    * Slope (S) Ultimately I want to create a regression that represents the cost function. Something like this:               fit <- lm(C ~ ß1xS + ß2xSD + ß3xTPA + ß4xVPT + ß5)


title: Autoclass in R/Python?

text: Are there any packages that implement the Autoclass/ Naive Bayes Clustering algorithm in R or Python? Alternatively, what are some other clustering algorithms that can handle both categorical and numeric variables that are implemented in either R or Python?


title: What type of analysis

text: I have two data sets to compare. Each is a list of billed amounts by diagnosis codes. The data differs in that the diagnosis codes may be different for some of the billed amounts. There are approximately 32 different diagnosis codes that were used. What type of statistical analysis is most appropriate in comparing these two datasets and why? I was told to complete a two way anova. Is that correct?


title: probability life similarities in a population cross section

text: I have just read up on Bouchard's Minnesota Twins study that turned up some amazing identical similarities in the lives of some identical twins reared apart. Both had same jobs, married and divorced women of same name, named dogs and children same etc. I'm not sure what percentage of their study group had these identical similarities(maybe 50% but not sure). I was just wondering if there have been any statistical studies comparing identical elements in lives of cross-sectional groups of two random people?


title: All neural network designs stop because of early stopping in MATLAB

text: I'm using `patternnet` for my binary classification in `MATLAB` and using `10-fold cross validation`. I'm repeating my neural network in a loop for some times to see the results. when i see neural network training window in `MATLAB` after 15~40 iterations the neural networks stops because of early stopping criteria (maximum 6 iterations is default value of early stopping in `MATLAB`). What do you think about this behavior? Is this normal or i should change something to improve it? I think in this condition the training can't train the network correctly, is that true? ( In all structures (inputs,neural network neurons number,layers,etc) i have this condition) Thanks.


title: Difference between Meta-Analysis, Meta-Regression and Moderator-Analysis

text: What is the difference between the classic meta-analysis (that aggregates effect sizes from a sample of studies to a summary effect size), meta- regression analysis and moderator analysis? As I understood moderator analysis is used to explain heterogeneity in a meta- analysis by regressing a vector containing the effect sizes of each study on specific variables that could explain the heterogeneity (e.g. publication year,...). But what's the difference between moderator-analysis and meta- regression? Thanks for your help!


title: Distribution for first time when the value is less than the previous one

text: Let $X_i, i \geq 1,$ be independent uniform (0, 1) random variables, and define $N$ by $$N=\min\{n:X_n < X_{n-1}\}$$. I need to prove that $$P\{N \geq k | X_0=x\} = \frac{(1-x)^{k-1}}{(k-1)!}$$. I am stuck when finding the quantities like $P(X_i \geq X_{i-1} | X_{j} \geq X_{j-1}, \forall j\lt i)$. I have also posted this in math exchange, but not convinced with the answer posted there. Another question is how to calculate $E[N]$ using conditioning on $X_1$.


title: Bayesian variable selection -- does it really work?

text: I thought I might toy with some Bayesian variable selection, following a nice blog post and the linked papers therein. I wrote a program in rjags (where I am quite a rookie) and fetched price data for Exxon Mobil, along with some things that are unlikely to explain its returns (e.g. palladium prices) and other things that should be highly correlated (like the SP500). Running `lm()`, we see that there strong evidence of an overparameterized model, but that palladium should definitely be excluded:               Call:     lm(formula = Exxon ~ 0 + SP + Palladium + Russell + OilETF +          EnergyStks, data = chkr)          Residuals:            Min         1Q     Median         3Q        Max      -1.663e-03 -4.419e-04  3.099e-05  3.991e-04  1.677e-03           Coefficients:                Estimate Std. Error t value Pr(>|t|)         SP          0.51913    0.19772   2.626 0.010588 *       Palladium   0.01620    0.03744   0.433 0.666469         Russell    -0.34577    0.09946  -3.476 0.000871 ***     OilETF     -0.17327    0.08285  -2.091 0.040082 *       EnergyStks  0.79219    0.11418   6.938 1.53e-09 ***      After converting to returns, I tried running a simple model like this                 model {         for (i in 1:n) {           mean[i]<-inprod(X[i,],beta)           y[i]~dnorm(mean[i],tau)         }         for (j in 1:p) {           indicator[j]~dbern(probindicator)           betaifincluded[j]~dnorm(0,taubeta)           beta[j] <- indicator[j]*betaifincluded[j]         }         tau~dgamma(1,0.01)         taubeta~dgamma(1,0.01)         probindicator~dbeta(2,8)       }      but I found that, pretty much regardless of the parameters to the chosen gamma distributions, I got pretty nonsensical answers, such as an unvarying 20% inclusion probability for each variable. I also got tiny, tiny regression coefficients, which I am willing to tolerate since this supposed to be a selection model, but that still seemed weird.                                             Mean        SD  Naive SE Time-series SE     SP         beta[1]       -4.484e-03   0.10999  0.003478       0.007273     Palladium  beta[2]        1.422e-02   0.16646  0.005264       0.011106     Russell    beta[3]       -2.406e-03   0.08440  0.002669       0.003236     OilETF     beta[4]       -4.539e-03   0.14706  0.004651       0.005430     EnergyStks beta[5]       -1.106e-03   0.07907  0.002500       0.002647     SP         indicator[1]   1.980e-01   0.39869  0.012608       0.014786     Palladium  indicator[2]   1.960e-01   0.39717  0.012560       0.014550     Russell    indicator[3]   1.830e-01   0.38686  0.012234       0.013398     OilETF     indicator[4]   1.930e-01   0.39485  0.012486       0.013229     EnergyStks indicator[5]   2.070e-01   0.40536  0.012819       0.014505                probindicator  1.952e-01   0.11981  0.003789       0.005625                tau            3.845e+03 632.18562 19.991465      19.991465                taubeta        1.119e+02 107.34143  3.394434       7.926577      Is Bayesian variable selection really that bad/sensitive? Or am I making some glaring error?


title: Using kriging with very sparse data

text: I have a very sparse set of measurements of a quantity, around 10 measurements at different locations across a continent. The distributions of locations is highly irregular. I would like to smooth these into values on a regularly spaced grid over the continent. We actually have a much denser set of observations of a related quantity and techniques are being developed to create a map of the first quantity using both sets of measurements. The intention is to estimate the accuracy of this hybrid approach using leave one out cross validation (testing each of the 10 sites in turn having created the map using all the others plus the additional measurements). As a baseline, I first need to see how accurate a map is, using LOOCV, created just from the original network of 10 stations. I should also point out this quantity varies over timescale of minutes and we want to produce a map in near real time. I am not an expert in geospatial techniques, however I comprehend the basics of Kriging and after checking the task view of CRAN, it seems autoKrige in the automap package would be a suitable function to employ. The results I am getting are nothing like what I expected so any help would be appreciated. Here is an example of one set of input training data (all available data for one particular time):                  lat   lon      value       -27.53 152.92   98       -35.32 149.00   79       -34.05 150.67   81       -12.45 130.95   92       -42.92 147.32   73       -22.25 114.08   91       -29.03 167.97  108       -31.94 115.95   89      As you can see, the target quantity has a range of values over a widely spaced area. I convert from a data frame containing lat, lon and the value into a spatialPointsDataFrame by doing               coordinates(xtrain) = ~lon + lat     coordinates(xtest) = ~lon + lat     proj4string(xtrain) <-CRS("")     proj4string(xtest) <-CRS("")      where xtest contains the grid of points I want to sample the map over. I then Krige this using               result <- autoKrige(value~1,xtrain,xtest)       I plot the results               plot(result)      and get the following image ![Kriging result](http://i.stack.imgur.com/VxXDX.png) As you can see, the result is essentially the mean value at all points. My understanding is that autoKrige should have determined the best values of the Kriging parameters and produced something more realistic than this. Just to illustrate, the range of the predicted values on the grid is 0.22, which is clearly much less than the 30 point scatter seen in the input values. So, where am I going wrong? Am I misusing the R package, or is my problem more fundamental in the way I am attempting to employ Kriging?


title: Come up with a simple experiment for linear regression

text: For my linear regression class we need to come up with some kind of simple experiment (we must collect the data, can't get it online) for which we can use the methods we have studied in class (simple linear regression, multiple linear regression, ANOVA, ...) to analyze the data. I am having the hardest time coming up with something. Anyone have a good suggestion? Last year someone did something silly like compare how fast mustard descends an incline compared to ketchup. EDIT: Thanks for all the helpful ideas guys :) I ended up with doing some physical experiments on myself on a stationary bike (max speeds/power output at different gear ratios). Probably not the greatest experiment, but it utilizes something I enjoy.


title: Unbalanced dataset - ROC curve to compare classifiers?

text: I use the machine learning software `WEKA` for data mining on biological data. I would describe my dataset as unbalanced: It comprises around 2000 instances, splitting in classes of 900, 500, 350, 160 that are very important to have in the dataset and some less important smaller classes that are nice to have but can be removed from the dataset if they confuse the learning to much.   Currently I am comparing many different classifiers. I am not a very experienced statistician, but I read that ROC curves are commonly used to evaluate the performance of machine learning classifiers. However, I also read that ROC has drawbacks when it comes to unbalanced datasets. Is there a better measure among the ones the `WEKA` output features (or can be calculated from them) for my dataset? Thats how the output looks like (here with the iris dataset):               === Stratified cross-validation ===           Correctly Classified Instances         144               96      %        Incorrectly Classified Instances         6                4      %        Kappa statistic                          0.94       Mean absolute error                      0.035      Root mean squared error                  0.1586     Relative absolute error                  7.8705 %     Root relative squared error             33.6353 %     Total Number of Instances              150                   === Detailed Accuracy By Class ===                         TP Rate   FP Rate   Precision   Recall  F-Measure   ROC Area  Class                      0.98      0          1         0.98      0.99       0.99     Iris-setosa                      0.94      0.03       0.94      0.94      0.94       0.952    Iris-versicolor                      0.96      0.03       0.941     0.96      0.95       0.961    Iris-virginica     Weighted Avg.    0.96      0.02       0.96      0.96      0.96       0.968               === Confusion Matrix ===             a  b  c   <-- classified as      49  1  0 |  a = Iris-setosa       0 47  3 |  b = Iris-versicolor       0  2 48 |  c = Iris-virginica


title: Generating correlated numbers from independent distributions

text: I'm trying to simulate this scenario: 10 different algorithms are solving a number of problems. All 10 are run on each problem instance, which means if a particular problem instance is hard, I expect all algorithms to perform quite bad and vice versa. Now, to simulate this, I created 10 different normal distributions from which I draw random numbers. These numbers represent the algorithms performance for some problem instance; however, since the distributions are independent, I'm not capturing the effect of the problem instance on the drawn numbers. Any idea how I can do this (i.e. correlated the random numbers that are drawn)? Regards,


title: How to sample from the product of a Gaussian and an inverse-Wishart distribution?

text: How do I sample from a distribution that is the product of a gaussian and an inverse-wishart distribution? I was going to use inverse transform sampling but a friend said that he thinks there is a much simpler way of doing it, a trick because the inverse-wishart is a conjugate prior of the gaussian. Does anyone know if there is an easy way to do this? Thanks very much!


title: How can I find correlations between crashes and system environments?

text: In my free time, I am working on a small web based system which collects crash reports (but not other, non-crashing bug reports) that are sent from Delphi Windows applications. For troubleshooting, users would love to have a data-mining feature to find relationships between hardware or operating system versions and the specific bug and/or crash. As an example how this should work:   * for every crash there is a report in the database, which has a fingerprint / hash code of the stack trace (call stack) at the moment of the crash to identify duplicates   * the algorithm checks if all duplicates of a bug report also have some other common attributes, for example a missing service pack of the operating system   * the analysis result lists all properties which the bug reports have in common Let's assume these automatic bug reports contain all key information like the names of all processes which are currently running, file names, version information of loaded DLLs, etc. How can I find correlations between repeated crashes and the environment? Are there specific algorithms or statistical methods that would help?


title: Analyzing ratios or numerator/denominator, with weights, measured over time at irregular intervals

text: I have gathered information on the rate increases allowed by a government regulator over about 17 years. As background, utilities in this area must apply to the regulator if they wish to raise their rates, in exchange for being granted a quasi-monopoly in the market. In this period, there are 65 'usable' decisions -- those with sufficient documentation. Each decision in a rate can be characterized by several bits of information:   1. The utility company    2. The utility company's industry (electricity, wastewater, natural gas, etc.)    3. The total revenue increase sought by the utility    4. The total revenue increase granted by the regulator    5. The % increase in revenues sought by the utility (implicitly, their total revenues)    6. The % increase in revenues granted by the regulator    7. The utility's rate of return on its operations in the status quo    8. The projected rate of return after the rate increase The 65 decisions are not evenly spaced over the 17-year period in question. To make matters potentially more complicated, some observations occur on the same date (four observations over two dates). From an exploratory point of view, I'm interested in detecting inflection points over time in these data, particularly in terms of the percent of the sought revenue increase ultimately granted by the regulator (#2/#1, above), the rate of return projected after the revenue increase, and the % increase in revenues granted by the regulator. I'm also interested in what (for lack of knowledge in this field) I'd call breakpoints or cleavages -- are there ways I could run a 2 sample t-test on all pairs of samples that have no chronological overlap? (e.g., sample 1 from time t0 to tn, sample 2 from t(n+1) to tf) From a research question perspective, I'm interested in whether the decisionmaking manifest in these data appear to have been impacted by a policy change that occurred at a known time, which (I would think), I could model as a categorical variable (0 = before change, 1 = after change) -- and if so, what was the sign of the impact, and was there a lag. Are there precautions I should be aware of with respect to the timed nature of these data? I consulted some materials on time series analysis but I'm a little confused as to what to do with irregularly spaced observations. As a final note, I guess I'd add that I'm working in R, but am flexible on that point.


title: Permutation test null and correct label distribution

text: Permutation tests can be used to test the likelihood of obtaining a single mean value (composed of the average of k-fold CV) classifier accuracy by chance, such as for a linear classifier using two classes. The permuted label distribution generates an empirical null distribution by randomizing the relationship between the labels and the instances. The exact p value is obtained by measuring how many values of the null distribution are equal or higher than the correctly labelled accuracy. My question is, what value should be ideally used as the 'correctly labelled accuracy'? For the next section, it is useful to refer to this illustrative plot. ![enter image description here](http://i.stack.imgur.com/XJIAc.png) For example, the null distribution is a spread around a mean value because because the random labelling of different conditions generates slightly different answers around chance on every repeat of the classification procedure (red distribution on plot). One repeat being the mean accuracy over k-fold crossvalidation. However, the 'correctly labelled accuracy' can also vary around a mean, simply because during k-fold cross-validation, different instances are assigned into different k-fold splits, resulting sometimes in good sometimes in better or worse values (green distribution on plot). This means that within this 'correctly labelled distribution' you will have some values that are extremely good, but highly unlikely (in the upper 97.5 confidence interval, see right hand P value). Yet often, studies report only a single accuracy instead of the mean of a distribution that has been repeated 10000 times (sometimes called bootstrapped mean and confidence intervals). So in the example plot, one could have as the mean 56 % and in the upper 97.5 tail 59 % correct, yet it is not clear where in this distribution the 'reported value' originates from. In summary: Question 1: If a permutation test is used to assess significance of a classifier accuracy, should the correctly labelled condition also be repeated also N- times to generate a correctly labelled distribution? Question 2: Where in the correctly labelled distribution should the compared value originate from?


title: Two-sample $t$-test vs Tukey's method

text: Even if the $F$ statistic in an ANOVA test is significant, we may need to do further testing before drawing conclusions. The most common method for doing this involves the use of a multiple comparisons procedure. I am trying to understand why Tukey's method is used over the two-sample $t$-test. I read the following below which suggested that doing all the pairwise tests (presumably it means two-sample $t$-tests) can lead to false positives (Type I errors). Why should this happen conceptually? > Tukey's procedure allows us to conduct separate tests to decide whether > $\mu_i = \mu_j$ for each pair of means in an ANOVA study of $k$ population > means. Like all multiple comparison procedures, Tukey's method is based on > the selection of a "family" significance leve, $\alpha$, that applies to the > entire collection of pairwise hypothesis tests. For example, when using the > Tukey procedure with a significance level of, say, 5%, we are assured that > there is at most a 5% chance of obtaining a false positive among the entire > set of pairwise tests, That is, there is at most a 5% chance of mistakenly > concluding that two population means differ when, in fact, they are equal. > This is very different from simply conducting all the pairwise tests as > individual tests, each at $\alpha=0.05$, which can result in a high > probability of finding false positives among the pairwise tests.


title: Conditions for the existence of a Fisher information matrix

text: Different textbooks cite different conditions for the existence of a Fisher information matrix. Several such conditions are listed below, each of which appears in some, but not all, of the definitions of "Fisher information matrix".   1. Is there a standard, minimal set of conditions?    2. Of the 5 conditions below, which can be done away with?    3. If one of the conditions can be done away with, why do you reckon it was included in the first place?    4. If one of the conditions can not be done away with, does it mean that those textbooks that did not specify it gave an erroneous, or at least an incomplete, definition?  >   1. _Zacks, The Theory of Statistical Inference (1971), p. 194._   >  The matrix $\mathcal{I}\left(\theta\right)$ is positive definite for all > $\theta\in\Theta$. >   2. _Schervish, Theory of Statistics (1997, corr. 2nd printing), Definition > 2.78, p. 111_   >  The set $C=\left\{x:f\left(x;\theta\right)>0\right\}$ is the same for all > $\theta$. >   3. _Borovkov, Mathematical Statistics (1998). p. 147_   >  $f\left(x;\theta\right)$ are **continuously** differentiable w.r.t. > $\theta_i$. >   4. _Borovkov, Mathematical Statistics (1998). p. 147_   >  $\mathcal{I}\left(\theta\right)$ is continuous and invertible. >   5. _Gourieroux & Monfort, Statistics and Econometric Models, Vol I (1995). > Definition (a), pp. 81-82_   >  $\frac{\partial^2}{\partial\theta_i\partial\theta_j}f\left(x;\theta\right)$ > exist > In comparison, here is the full list of conditions in _Lehman & Cassella. Theory of Point Estimation (1998). p. 124_: >   1. $\Theta$ is an open interval (finite, infinite, or semi-infinite) >   2. The set $C=\left\{x:f\left(x,\theta\right)>0\right\}$ is the same for > all $\theta\in\Theta$. >   3. $\frac{\partial f\left(x;\theta\right)}{\partial\theta_i}$ exists and > is finite. > And here is the complete list of conditions in _Barra, Notions fondamentales de statistique mathematique (1971). Definition 1, p. 35_ : > The score is defined for all $\theta\in\Theta$, each of its components is > square-integrable and has integral $=0$. It is interesting to note that neither Lehman & Cassella nor Barra stipulate that $\int f\left(x;\theta\right)\space \mu\left(dx\right)$ be differentiable under the integral sign w.r.t. each $\theta_i$, a condition that occurs in most other textbooks i surveyed.


title: Bayesian model with unknown mean and variance with lognormal prior

text: For $i=1, \ldots, K$ and $j=1, \ldots,n$, assume the following model. \begin{align} X_{ij} \mid \mu_i, \sigma^2 & \stackrel{_\text{iid}}{\sim} N(\mu_i, \sigma^2) \nonumber \\ \mu_i & \stackrel{_\text{iid}}{\sim} N(\mu, \tau^2) \nonumber \\ \ln(\sigma^2) & \sim N(\mu_v, \tau_v^2) \end{align} First assume all the prior parameters are known. Since this is not the conjugate form, is there a way to estimate the posterior $\pi(\mu_i \mid \mathbf{X_i})$ using a closed form density? Thank you very much.


title: Summarize contingency table in R

text: I know table can count contingency table, what if I want to summarize all elements that fall in a cell of the table using some type of function, such as mean? For example               facA = rep(1:3, 4)          facB = rep(c(rep(1, 3), rep(2, 3)), 2)          value = rnorm(12)      How do I summarize the mean of 'value's corresponding to facA:facB and make it as a table?


title: huge sample sizes and null hypothesis testing

text: I'm kinda brand new to statistics and I'm developing classification algorythm. My method is based on simple chi-square goodness of fit. I am counting the effect size of known cases to predict the future ones. My problem is that I don't really know how to deal with large sample sizes, neither do I know when is a sample size "too big". Did research the answer, but did not find solution. So what is the best approach to take when you have a double possible outcome and like 2000-3000 known cases where you know the outcome? Obviously it isn't safe for NHST as large sample sizes will yield false positives (as much as I know). Can anyone suggest some good approach to take, or maybe a good article on the topic? Thanks for help!


title: How can I produce a plot showing the directional angles of my points?

text: I want to be able to compare the angles of neighbours in a herd of hippos. I have data for the x and y coordinates and the angles that they are facing (using imageJ, angles are between -180 and 180 with 0 being a horizontal straight line across the middle of the image). I have plotted the $g(r)$ function to see whether my hippos 'like' to orient themselves towards individuals that are nearby, but I would love to be able **to compute a visual representation of their orientation in the herd** using `R`. Any help would be much appreciated, I'm not really `R` literate yet!


title: How to sample natural numbers, such that the sum is equal to a constant?

text: Say I have $N$ items that are partitioned / clustered and I want to randomly repartition these items, such that the distribution of sizes of the clusters is 'similar' to those that I already have. I'm viewing this (perhaps unhelpfully), as trying to sample $k$ natural numbers, such that $x_i \geq 1$ and $\sum_{i=1}^{k}x_i=N$. $N$ in this context would probably be in the thousands and the distribution of the current cluster sizes would be quite skewed, with a small number of large clusters and a large number of small clusters. I've looked at hypergeometric, binomial, geometric distributions, but none of these seem to quite fit what I'm looking for - I'm guessing that it's more complicated than a simple distribution and some type of Markov-like process would be needed. Anybody have any ideas?


Let $X_{i}, i\ge 1$, be i.i.d. random variables defined on a probability space $(\Omega, \mathcal F,P)$ such that $P(X_{i}=1)=P(X_{i}=-1)=\frac{1}{2}$. Consider the filtration $\mathcal F_{n}=\sigma(X_{1},\dots,X_{n})$ on this space and the random walk $R_{i}=\sum_{i=1}^{n} X_{i}$. Show that $R_{n}^{2}-n$ and $(-1)^{n} \cos(\pi R_{n}) $ are $\mathcal F_{n}$-martingales.


title: Ranking topics in K-Means

text: I understand that clustering is meant to group items together, is there any ways that we can quantify saying Cluster A is more important than Cluster B? Other than counting the number of items in a cluster?


title: self organizig map with spectral clustering

text: I'm looking for implementation of **self organizing map** which for the clustering part **uses spectral clustering rather than k-means**. Does anyone knows something about it ?


title: Writing R Function To Determine Which Statistical Test Should Be Run

text: I am trying to write a function in R that determines which of the following statistical tests should be run on a set of data: {paired-t, sign, two-sample t, pooled two-sample t, two proportions, F} -Data sets have two columns, representing the X and Y values. -The focus of the test on the data sets are one of the following types: {center, spread, count} (for count the data in columns 1 and 2 will be 1’s and 0’s for the number of successes and failures respectively). -Data sets are marked as either dependent or independent. I am comfortable with coding in R and I have written the functions to determine if the data set is normal and if there are outliers, but I am not clear on how, given this information I could determine which of these tests should be run. Could someone please explain how these criteria suggest which of these tests should be run?


title: two sample mean comparison with small and nonnormal data

text: I have three groups. Both have small sample sizes (<10). And we can't check the assumption of normality. I want to compare any two of them. What to do in this case( unequal small sizes, non-normality)


title: Kolmogorov-Smirnov test with estimated parameters

text: I've read through some past posts about using the KS test when the distribution parameters are estimated. This one in particular was very helpful, and I implemented Greg Snow's suggestion. However, I am not clear on what to do when testing my proposed distribution against a second set of data. This is my procedure:   1. Collect a sample of data over a given time period   2. Determine the parameters of a hypothesized distribution from the data, and test the goodness of the fit with the KS test (modified as per Greg Snow's answer in the other post)   3. If it passes the KS test at an acceptable confidence level, test the distribution (with the calculated parameters) against another set of data collected at a later time. In step (3), I'm interested in how well the _particular_ parameters values calculated in (2) work at predicting future data. That is, I'm now treating the previously calculated parameter estimates as fixed inputs. Hence, it seems to me I should be using the standard KS test with the canonical distribution for the test statistic. Is this correct?


title: why i do not have sgplot with runing this program in sas 9.2

text: data a;     input treat$    rep boreria mimiosa asysta  axono   seleria paspa   macran  melas   scopo   cleom   otto;     cards;     M   1   0.20    0.00    0.00    0.00    0.00    12.20   0.00    0.00    0.00    4.90    0.00     M   2   0.00    0.00    0.00    0.64    0.00    0.00    0.00    0.00    0.00    0.00    0.00     M   3   0.00    0.00    0.00    0.23    0.00    0.00    0.32    0.00    0.00    0.00    0.00     A   1   0.00    0.00    7.80    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00     A   2   0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00     A   3   0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00     CA  1   9.40    1.24    1.59    5.20    0.00    1.77    0.00    0.00    0.12    0.00    0.00     CA  2   0.75    0.00    0.00    2.10    0.00    0.00    0.00    0.00    0.00    0.00    2.66     CA  3   2.00    0.99    0.00    1.10    0.00    0.00    0.00    0.76    0.00    0.00    0.00     P   1   0.00    0.00    0.00    1.56    0.97    0.00    0.00    0.10    0.00    0.00    0.00     P   2   0.02    0.00    0.00    0.93    0.00    3.50    0.00    0.00    0.00    0.00    0.00     P   3   0.00    0.00    0.00    0.00    0.00    6.20    1.00    0.00    0.00    0.00    0.00     W   1   13.40   15.80   0.30    104.50  0.00    127.70  0.00    0.00    0.00    0.00    0.00     W   2   4.80    0.00    0.64    121.10  41.10   58.10   0.00    9.20    0.00    0.00    0.59     W   3   0.00    0.00    0.00    53.00   0.00    54.70   0.00    10.10   0.00    0.00    0.00     ;     data b;     set a;     array log{*} boreria--otto;     do i = 1 to dim(log);     log(i)=log10(log(i)+1);     end;     output;     drop i;     run;     ods output  Eigenvectors=vector ;     proc princomp out=pcacomp  ;         VAR boreria--otto;     RUN;          data vec;     set vector (rename = (variable=Weed prin1=vec1 prin2=vec2));     vec1 = vec1*7;     vec2 = vec2*7;     drop prin3-prin10;     run;          data new;     merge pcacomp vec;     run;                    ods listing style=statistical image_dpi=300;     proc sgplot data=new;     VECTOR x=vec1 y=vec2/group=weed ARROWHEADSHAPE=filled DATALABEL=weed;     scatter x=prin1 y=prin2 / group=treat name='site';     ellipse x=prin1 y=prin2 / LINEATTRS=(pattern=2) TRANSPARENCY=0.1 ;      refline 0 /axis=x LINEATTRS=(PATTERN=1);     refline 0 /axis=y LINEATTRS=(PATTERN=1);     xaxis label='PCA1 (32.0%)' grid ;     yaxis label='PCA2 (18.2%)' grid;     KEYLEGEND 'site' /title='Treatment' LOCATION= inside position= topLEFT across=3 noborder;      run;


title: Monte Carlo / Nonparametric confidence intervals for mean estimate

text: I have 2x2 design with n = 3 (averages) for each group (see here). I am not sure that my data follows any particular distribution, but I would like to use confidence intervals as error bars in a dot plot. Thus I can't really use sd, se or bootstrapping from a distribution. The approach I am after is described in this article. My initial feeling is that I should use boot package. I have 3 replicates, which consist of different number of measurements. I would like to use averages of these replicates to calculate the confidence intervals. Something like:               x1 <- rnorm(50, 14,3)     x2 <- rnorm(35, 7,1)     x3 <- rnorm(40, 15,9)          d <- c(mean(x1), mean(x2), mean(x3))     b <- boot(d, function(u,i) mean(u[i]), R = 999, sim = "ordinary")     boot.ci(b, type = c("all"))      The help file for `boot()`function says "For the nonparametric bootstrap, possible resampling methods are the ordinary bootstrap, the balanced bootstrap, antithetic resampling, and permutation". Without background in statistics, I am having hard time understanding what do these different options do and what would be the most suitable for this case. `boot.ci()` function gives 5 different answers, which all seem to differ considerably. I am pretty confused. Could anyone clarify following related questions:   1. Is this a right approach?   2. Which combination should I use in this case?   3. Why so many options?


title: Distribution of three dependent random variables

text: I have three dependent random variables $X_1$, $X_2$, $X_3$, ($X_i\in\left[0, 1\right]$) that satisfy additional condition $\sum_i{X_i} = 1$. If (marginal) distributions of $X_i$ ($i = 1, 2, 3$) are given, I need to determine likelihood of a certain combination of those variables. More formally, if $a_1, a_2, a_3\in \left[0, 1\right]$ such that $\sum_i{a_i} = 1$ and $\epsilon > 0$, I would like to know the probability: $$ P((\forall i) |X_i - a_i| < \epsilon). $$


title: What happens when SEM is conducted with a rank-based DV?

text: Can SEM be conducted with a rank-based dependent variable? If so, which software package would be the way to go? Obviously, you'll get LISREL output and such, but I'm concerned as to how the rank-based DV would affect parameter estimates and was wondering if something can be done to combat any bad things that come from ignoring the distributional assumption of the DV. Thanks!


title: Generating samples from Copula in R

text: Suppose I want to model dependence between $d$ r.v.´s $Y_1,...,Y_d$ with the copula $C_\theta$, where $\theta$ are the corresponding parameters of that copula. I've also determined the correlation using Kendall's tau $\tau(Y_i,Y_j)=\gamma$ for all $i\neq j$. Now I want to generate samples from this joint distribution: Am I on the right track when I say that, from Sklar's Theorem, we have that if $Y_k$'s corresponding distribution function is $G_k(y_k)$ and we have $X_1,...,X_d$ with distribution functions $F_1(x_1),...,F_d(x_d)$. Then $\bf{Y}$$=\left(G_1^{-1}(F_1(X_1)),...,G_d^{-1}(F_d(X_d))\right)$. So say for example that I want to generate samples of our t-distributed random vector $\bf{Y}$ from the join distribution using a Gaussian copula, I could just simply generate the normal random variables $Y_1,...,Y_d$ (independent) and then apply the t distribution function on the r.v. and then take the normal quantile function of that, using the $\Sigma$ where $$\Sigma_{i,j}=\begin{cases} 1, & \text{if $i=j$}.\\ sin(\frac{\pi}{2}\gamma), & \text{otherwise}. \end{cases}$$ as one of the parameters. Or should the dependence be taken into account for when I model the t-distributed vector $\bf{Y}$? I guess the first makes more sense since the copula is supposed to "model the dependence". Say instead I want to do the same using the copula package in R. Is it correct to first create the copula object with $\Sigma$ (as before) using _normalCopula_ and then _mvdc_ and _rmvdc_ to generate the samples? But what are the parameters for the mvdc function? The example code from R's web page gives the following: mvdc(normalCopula(0.75), c("norm", "exp"),list(list(mean = 0, sd =2), list(rate = 2))) But shouldn't the argument in the normalCopula function be a 2x2 matrix since its bivariate. Also, say now that I have $d$ equally t-distributed r.v.'s and I dont want to write c("t", .... ,"t") etc. How can I write this neatly? As you can see I'm slightly confused about the parameters in the copula and in what step the correlation should be used (I'm also not sure when it is eligible to say that $\rho = sin(\frac{\pi}{2}\tau)$), so if someone could help clarifying that I'd be very thankful. Also, there are lots of examples of 2-dimensional (and some on 3-dimensional) examples of this on the internet since they are more illustrative than for higher dimensions. So if someone have references to examples of modelling with copulas in higher dimensions in R I would gladly accept that. * * * Edit: Ok, say I want to model the joint distribution of two financial assets where the marginals are t-distributed with 3 degrees of freedom and standard deviation 0.05. Their correlation is estimated with $\tau(Y_1,Y_2)=0.5$ I've assumed their joint distribution should be normal, so I use a normal copula. To generate samples from their joint distribution. I first generate samples from the normal distribution: X<-rmvnorm(N,mu,S) #mu is zero vector and #S is as $\Sigma$ before with $\gamma = 0.5$. ![Sample from Multivariate Normal Distr](http://i.stack.imgur.com/WrMgp.jpg) Now I transform each sample point $(X_k^1,X_k^2)$: Y<-qt(pnorm(X),df=3) ![enter image description here](http://i.stack.imgur.com/ZgoR8.jpg) Is this correct? But how do I control the variance/standard deviation of the marginals? Do I do that when I generate the multivariate normal or in the qt function?


title: Clustering text with python

text: I asked on stackoverflow but they suggest me to move here for better answers. I copy paste the question. I decide to play a little with similarities and clustering text. I have already create the tf-idf and symmatrix matrix of similarities. Now I want to implement something for clustering into groups. I made a research and I have found hcluster (Hierarchical clustering) and k-means libraries. Which of those is better in base of accuracy? Do you know any better way even if there isn't a ready library. I could write the code if I knew the algorithm. Also, this way is O(n^3). If I would like to sacrifize a little of accuracy to win on calculating time, do you have something to suggest?


title: How exactly does the R random number generator work?

text: When I generate a random number using the following procedure:               set.seed(123)     rnorm(1)      I would like to know what exactly is happening in this case? Is R reading from a table of random numbers? What exactly is 123 referring to? Is it possible to access this table to see how it looks?


title: How to interpret a non significant independent variable?

text: I conducted a regression analysis using R's lm() function. One of the independent variables shows no significance (p = 0.89), which contradicts the hypothesis that is should have a significantly positive effect on the dependent variable. How do you interpret that? Can you say that it has no positive effect on the dependent variable, just because it is not significant - even though it is not significantly negative?


title: Distinguishing statistical (global) and network (local) effects

text: I am analyzing how users of specific service affect each other by observing communication between them and changes in membership plans. The social network consists of 3M users and 40M connections among them. To determine the influence of a specific user I observe the sequence of membership plans changes. For example there is a network of users A, B, C, D, and E and the connections               A SENT_MESSAGE_TO B     C SENT_MESSAGE_TO B     D SENT_MESSAGE_TO B     D SENT_MESSAGE_TO E      and the membership log               B changed the subscription plan from free to premium in 2011-10-11     A changed the subscription plan from free to premium in 2011-10-18     D changed the subscription plan from free to premium in 2011-10-22      As the users are connected to each other there are (I assume) two possible causes for membership plan change:   * The person was either influenced by advertisment, personal satisfaction or in generally the global effects had a major influence on him    * the person noticed the behavior of his friends (connections) and imitated the majority or a specific neighbor As there is possible to observe the subscription plan changes in whole network I guess I could statistically determine, what was the effect of a particular advertising campaign or price change (for example I know the dates of each campaign and price change) in global scale but don't know how to use this insight to distinguish interpersonal effects from global effects. One guess (a simple example) would be to say 'if the time difference in plan change between two connected users is more than [determined period] then we assume the second change wasn't affected by first one'. Another guess could be 'if two users exchanged less than [specific number of messages] then there is no influence among each other'. Anyway I would like to find or develop a solution for distinguishing global and network effects in social network data and would be very thankful if anyone suggested any good article with similar problem or suggest an approach for dealing with such problems. Thank you!


title: Negative binomial -- IRR interpretation for predictors

text: I have a zero-inflated negative binomial model. I have used incidence rate ratios and I'm trying to interpret the coefficients in relation to my predictors. Most of my predictors are continuous variables of census data -- ie: % of the population that is Hispanic; % of the population less than age 18, etc. I know that the IRR is normally interpreted as the rate ratio for a 1-unit increase in the independent variable, but what does this mean in terms of these continuous predictors -- does this mean the IRR is the estimated rate ratio for a 1% increase in % Hispanic. Is there a way I can scale this so it can be interpreted to be the estimated rate ratio for a 10% increase in the % Hispanic? Also, one of my IRR's is 20. Does that seem unusually high?


title: Compute p-value in paired bootstrap

text: I ran across a new paper from the Berkeley NLP group on statistical testing, An Empirical Investigation of Statistical Significance in NLP. There is pseudocode for computing a p-value in the paper, basically, the idea is that the sample set of $x_1,x_2,...,x_N$ are sampled with replacement from data $x$. Then $\text{p-value} = \text{count}(\delta(x_i) > 2\delta(x))/N$, where $\delta(x_i)$ is a metric gain. I could understand the formula for computing the p-value in Koehn's paper Statistical significance tests for machine translation evaluation, in which: $\text{p-value} = \text{count}(\delta_a(x_i) < \delta_b(x_i))/N$, where $\delta_a$ and $\delta_b$ are the metric gain for system $a$ and $b$ respectively. Is there any explanation or reference for the formula $\text{p-value} = \text{count}(\delta(x_i) > 2\delta(x))/N$. The authors also noted that if the mean of $\delta(x_i)$ is $\delta(x)$ and $\delta(x_i)$ is symmetric, then both formulas above are equivalent.


title: How to visualize counts of different categories

text: Suppose we have $10$ boxes and we are interested in measuring the number of apples, oranges and pears in each. What is a good way to visualize how the boxes relate to each other in terms of the distribution of fruits? For example, box 1 might have 5 apples, 10 oranges, and 6 pears while box 10 might have 1 apple, 1 orange and 1 pear. **Added.** Is looking at a distance matrix appropriate? Is there a way to visualize a distance matrix?


title: What test (or algorithm) should I use for the following use case?

text: I have a distribution of visitors to a website on Monday to Sunday like this:               M     T   W  Th   F   S  Su     345 467 560 350 430 689 490      I know that if I have a hypothesis of distribution that I suspect such as (M 20%, T 30%, W 10%, Th 5%, F 5% S 20% Su 10%) etc, I can use a Chi Square distribution to figure out if my suspicion is correct - as explained in this video, Pearson's Chi Square Test (Goodness of Fit), on Khan Academy. My question is how do I figure out an acceptable distribution given an observed data. I think since I am not after a True or False, rather a set of numbers, I suspect I am after some kind of algorithm that gives out a range of percentages that specify a region which are all acceptable under a chi-square test. Any pointers greatly appreciated!


title: Inferring from a combination of uncertain and certain data

text: I am trying to estimate the surface (isochrone), $z_i(x)$ for which $T(x,z)=0$ from noisy measurements of $T(x,z)$ everywhere and 3 almost noise free control points: ![isochrone](http://i.stack.imgur.com/6eDJG.png) Instead of using the $T(x,z)$ data directly I first estimate the gradient: $T_x$ and $T_z$ by finite differences on $T(x,z)$. I now have two pde's to solve: $ \frac{\partial}{\partial x}T=T_x, \frac{\partial}{\partial z}T=T_z$ Using finite differences again I express the derivation operators as matrices: $\ D_x T=T_x, D_z T=T_z $ where $T_x$ and $T_z$ are now column vectors made by concatenating all the columns of their respective matrices vertically, and $D_x$ and $D_z$ are matrices that when multiplied with the $T$ vector gives the partial derivative with respect to x and z. In addition to these pde's I also have equations for the 3 control points: $ T(xc_i,zc_i)=0$ (Weighted) Least Squares Error (LSE) is a framework that allows me to mix certain and uncertain data. My (weighted) linear equations are now: $\ D_x T=T_x $ $\ D_z T=T_z $ $\ \alpha T(xc1,zc1) = 0$ $\ \alpha T(xc2,zc2) = 0$ $\ \alpha T(xc3,zc3) = 0$ where I have introduced $\ \alpha $: a parameter that decides how much the control points are honored relative to the measurements of $T(x,z)$. **Is there some way to estimate $\alpha$ or alternatively express it as a more intuitive expression?** So far I have just tuned this parameter and compared the offsets from the control points iteratively. It would be good to have some other way to set this, e.g. max allowed z offset from control point to resulting isochrone. **Do statisticians have some other way than LSE of combining certain and uncertain data?** Is Bayes such a method? If so what is the a priori in my example?


title: Comparing the effect of a treatment that was optional for its receivers

text: In 2012 we collected data at our university about retention of students from first semester to second semester along with some other variables. The retention variable is binary, 'retained' or 'did not retain'. In 2013 we introduced a new system. Student success advisers were assigned and any student when needed could go to them. We collected some new variables associated with the introduction of the student success advisers together with those collected in 2012. Now, going to these advisers was not mandatory. Those who felt a problem could go to them. Just looking at 2013 data alone makes it seem that the advisers make things worse: the students they see do worse, and are more likely to leave university. Only the students who faced problems probably went to the advisers and those who were doing well probably didn't. If 10 students went to the advisers with problems and 8 retained to semester 2, that is surely a success. But if 80 out of 100 students who didn't have any problem and thus didn't go to the adviser retained to semester 2, that will have log odds similar to those who went to the advisers and thus tell you that going to the advisers had no significant effect! Part of the reason was a major policy change that affected how universities in our country took in students. There was a "cap" in place before 2013, a limit on how many students a university could take on and expect government funding for. In 2013 that was removed. So, successful, desirable universities started opening their doors to more students, and more students were able, suddenly, to get into desirable universities. Our university is not so desirable. So what probably happened was that our university was forced in 2013 to take significantly "worse" students, which could well lead to greater retention difficulty. What my bosses think is that we need to weight the 2012 data so that it matches the 2013 data on some of the variables that may have been sensitive to the removal of the cap on student entry (or weight the 2013 data so that it matches the 2012 data). Then we can compare two different cohorts as if they were a proper control group for each other. In fact, this is illusory. Some suggested using SPSSINC RAKE, but I am not getting any clue why. What we want to find out (very broadly) is whether introducing the advisers etc. worked, and the degree to which it worked.


Let $X$ be an integrable random variable defined on probability space $(\Omega , \mathcal F,P)$, and let $\mathcal F_{n},n\ge0$ , be a filtration on this space.   show that $X_{n}= E[X|\mathcal F_{n}]$ is an $\mathcal F_{n}$-martingale.


title: Differences Between Logistic Regression in Statistics and in Machine Learning

text: I just found out that machine learning also has logistic regression as one of its methods. Can someone please tell me the differences between logistic regression in statistics and machine learning? I've seen lecture slides on logistic regression from a machine learning course, but I can't see the difference with the coverage of logistic regression in a statistics course. Does logistic regression in machine learning have no need to check for multicollinearity? The reason I asked this is because I've tried to run a dataset through R's `glm` function with binomial logit, and then I ran the same dataset through Apache Mahout's `trainlogistic`. But the resulting coefficients are different. This is the command I use in R:               w1.glm <- glm(anw ~ cs, data = w1, family = "binomial")      This is the result of `summary(w1.glm)`:               glm(formula = anw ~ cs, family = "binomial", data = w1)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.5400   0.1073   0.1924   1.0047   1.0047            Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  0.42077    0.02588   16.26   <2e-16 ***     cs           1.89342    0.06427   29.46   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 11762.5  on 10660  degrees of freedom     Residual deviance:  9250.3  on 10659  degrees of freedom      And this is the command I use in Mahout:               /usr/local/mahout/bin/mahout trainlogistic --input w1.csv --output ./model --target anw --categories 2 --predictors cs --types numeric --features 20 --passes 100 --rate 50          Running on hadoop, using /usr/local/hadoop/bin/hadoop and HADOOP_CONF_DIR=     MAHOUT-JOB: /usr/local/mahout/mahout-examples-0.8-job.jar     20     anw ~      -19.553*cs + -7.512*Intercept Term                 cs -19.55265           Intercept Term -7.51155         0.000000000     0.000000000     0.000000000     0.000000000     0.000000000     0.000000000   -19.552646543     0.000000000     0.000000000     0.000000000     0.000000000     0.000000000     0.000000000    -7.511546797     0.000000000     0.000000000     0.000000000     0.000000000     0.000000000     0.000000000      13/11/01 02:04:47 INFO driver.MahoutDriver: Program took 22118 ms (Minutes: 0.3686333333333333)      **Edited:** Added the reason I asked the question in the title. Added the commands used to execute `glm` in R and `trainlogistic` in Mahout.


title: Caret with rxDForest as custom model

text: I'm trying to use the `train` function from the caret package to tune the parameters of the `rxDForest` from RevoScaleR package (I can't use the random forest because my dataset is too big), but I can't manage to make my code work with resampling:               library(caret)     library(doParallel)     library(pROC)          registerDoParallel(cores=4)          # ----------------     # Import data train          load(file = "data_model_rf.RData")          # ----------------     # tune rf          rxDForest <- list(type = "Classification", library = "RevoScaleR", loop =     NULL,                       parameters = data.frame(parameter = c("mTry",     "maxDepth"),                                                    class = rep("numeric", 2),                                                   label = c("mTry", "maxDepth")),                       grid = function(x, y, len = NULL) {                                       p =     floor(sqrt(ncol(x)))                                       expand.grid(mTry =     c(p/2,p,p*2),                                                    maxDepth = c(1,3,5))},                       fit = function(x, y, wts, param, lev, last,     weights, classProbs, ...) {                                      tmp <-     as.data.frame(x)                                      tmp$Class <- y                  f <-     as.formula(paste0("Class", " ~ ",                       paste(colnames(tmp)[-c(ncol(data_model))],collapse = "+")))                  rxDForest(f,     tmp, maxDepth = param$maxDepth,                                      mTry = param$mTry,...)},               predict = function(modelFit, newdata, preProc     = NULL, submodels = NULL) {                                     rxPredict(modelFit, newdata)},               prob = function(modelFit, newdata, preProc =     NULL, submodels = NULL) {                   out <-     rxPredict(modelFit, newdata, type="prob")                   colnames(out)     <-  modelFit$obsLevels                                           out[,1:length(modelFit$obsLevels)]},               sort = function(x) {x[order(x$mTry, x$maxDepth     ),]},               varImp = function(object, numTrees = NULL,     ...) {                        varImp <-     object$importance                  out <-     data.frame(varImp)                      colnames(out) <- "Overall"                      rownames(out) <- rownames(object$importance)                                         out})          rfGrid <- expand.grid(mTry = c(14), maxDepth = c(5))     fitControl = trainControl( method = "repeatedcv", number = 10, classProbs =     TRUE,                                 repeats = 3, summaryFunction = twoClassSummary)               # ----------------     # tune rf          rfFit <- train( data_model[,1:(ncol(data_model)-1)],     data_model[,ncol(data_model)],                     method = rxDForest, trControl = fitControl, tuneGrid =     rfGrid,                     metric = "ROC", verbose = TRUE, nTree = 20,     importance=TRUE )      I got the following error:               Error in train.default(data_model[, 1:(ncol(data_model) - 1)], data_model[,  :      final tuning parameters could not be determined     4: In nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,  :     There were missing values in resampled performance measures.     5: In train.default(data_model[, 1:(ncol(data_model) - 1)], data_model[,  :     missing values found in aggregated results      The code works if I switch to `method = "none"` in the `trainControl` function, so it's not a problem with the `fit` or `predict` functions. Thanks for your help. EDIT : Documentation for the rxDForest : https://docs.google.com/file/d/0B-MIQzwMSpnxSlZXNjIxektEWmc/edit?pli=1 rxDForest accept only a specific "rxFormula" which doesn't include class~. but works with my syntax. I've tried the new grid matrix (included x = matrix(1, ncol = 9, nrow = 10) in the grid function) but it doesn't work. I'm using the tunegrid option so the grid function is not supposed to be used, isn't it ? Thanks.


title: Problem with $\chi^2$ test

text: I am very new to statistic analysis and I am having some problems with my data. Basically, my categorical variables are the following: grammatical structures (four of them) and languages (five of them). The total number of these structures is 1218. I want to answer the research question if the observed preference shown by some languages for certain grammatical structures is merely due to random chance or there is a relationship between them. I applied the chi-square test and I obtained the following p-value: 1.5168E-140 (chi-square value = 692.813, due to a huge difference between observed and expected frequencies, and a df = 12). I am really puzzled by this result. My guess is that the data are too dispersed and the chi-square test is inconsistent here. So does that actually mean that the sample is biased or that it does not make much sense to answer that particular research question in the beginning? Or maybe should I use some other test? Thanks in advance.


title: Time line analysis

text: I am conducting research on the relationship between a person's birth order and later risk of obesity using data from several 1-year birth cohorts (e.g. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2908417/). A key challenge is that birth order is linked to other features such as maternal age, number of younger and/or older siblings, and birth spacing, that may also influence the outcome via different mechanisms. Further, any influence these things have on later obesity risk might be modified by the gender composition of the siblings, including the "index child" (the participant in the birth cohort). For each index child, one could draw a timeline that showed all of the births in the family, with maternal age at the time variable. ![enter image description here](http://i.stack.imgur.com/s90mj.png) I am trying to identify methods to analyze these kinds of data, where the order, timing, and nature of events might all be important. I am asking this question here because of the diversity of applications members work with - I expect someone has some immediate suggestions that will take me much longer to identify alone. Any nudges in the right direction(s) would be greatly appreciated. Related question(s): How should I analyze data on birth intervals of women?


title: Grouping observations based on variables that sum to one

text: I have a problem where I am trying to group observations (most likely using k-means or a similar unsupervised learning tool) where each observation includes _n_ -variables, with the total sum of these variables equal to one. We are therefore grouping observations, based on the probabilities of each of _n_ potential states or outcomes for each observation. **For example** : if we are testing the manufacturing of dice, and we take a sample of 1000 dice, each die would be one observation and we would record V1 = the percentage of the time that a one appeared when we rolled that particular die, V2 = % two, etc. We would then cluster the die based on the percent of the time that each die appeared as a 1, 2, 3, etc (to see if, for example, one of the die machines was improperly). While we could look at the percentages, it is perceived that clustering these observations will help point to any under-lying trends (colors, materials, machines, etc) that would not be apparent without the unsupervised learning techniques. I know this is a bit of a strange problem, but I'm sure some literature exists for it. What is this type of problem called? Can you point me in the direction of some work that has been done to study this type of problem??


title: Mixture of binomial distributions

text: I am experimenting with a mixture of binomial models. Consider a binary variable $y_i$. Furthermore, there are two sub-groups in the population (not known a priori and not observable): $z_i=0$ or $z_i=1$. For one sub-group of the population, the following holds true: \begin{equation} P(y_i=1|z_i=0)= \xi, \end{equation} with $\xi$ a fixed probability (not estimated by the model but provided by the user). For the other sub-group, we have: \begin{equation} P(y_i=1|z_i=1)= g_i = g(x_i' \theta_g) = \frac{1}{1+\exp(-x_i'\theta_g)}, \end{equation} with $x_i$ a vector of independent variables (including intercept), $\theta_g$ a parameter vector and $g_i$ short-hand notation. The probability that an observation belongs to a certain sub-population is equal to: \begin{equation} P(z_i=1)= h_i=h(\tilde{x}_i' \theta_h) = \frac{1}{1+\exp(-\tilde{x}_i'\theta_h)}, \end{equation} with $\tilde{x}_i$ a vector of independent variables and $\theta_h$ the corresponding parameter vector. Note that, in general, $\tilde{x}_i$ may be different from $x_i$. In a special case, $\tilde{x}_i$ only contains an intercept, in which we have fixed mixing proportions (I have experimented with both, and both lead to the same problem, described below - let variant 1 be the one with intercept only, and variant 2 the one with more independent variables). The problem occurs when trying to estimate this model. I have implemented the EM algorithm (double/triple/quadruple checked, I am 95% sure the math is correct), but it seems my model is not identifiable (neither of both variants). When creating a contour plot (or plotting the surface) for two independent variables, there is no 'peak' but rather a 'ridge' (see below image). ![Contour plot](https://dl.dropbox.com/u/46839575/contourplot.JPG) I can intuitively understand this when looking at: \begin{equation} P(y_i=1) = h_i\cdot g_i+(1-h_i) \cdot \xi. \end{equation} I can imagine that one could increase $h_i$ and decrease $g_i$ in the right proportions and arrive at the same probability $P(y_i=1)$. And then probably this holds true for the LL as well. I have started to look in the literature for solutions for this problem, but I did not find a final answer yet. Is there a solution for this identifiability problem (at least I think that is the problem)? Could some additional constraints solve the problem, or is there another solution? I it, after all, possible to create a mixture of binomial distributions? I hope my explanation is clear. Thanks in advance for your help. Best regards, Thomas


title: When an ARMAX model is stationary? Why stationarity or invertibility is needed?

text: Let $y_t$ a stochastic process and $\tau_t$ presents the time duration between the $t$ and $t-1$ event.The ARMA(p,q,r) with exogenous variables is defined as: $$ y_t = \varepsilon_t + \sum_{i=1}^p \alpha_i y_{t-i} + \sum_{i=1}^q \theta_i \varepsilon_{t-i} + \sum_{i=0}^b \eta_i \tau_{t-i}.\,$$ where $\eta_1, \ldots, \eta_b$ are the parameters of the exogenous input $\tau_t$ and $\varepsilon_t$ is a white noise. # Question 1: I have read that this process is stationary if and only if the roots of the AR polynomial must be outside the unit circle, but I don't understand it. I understand the proof when the model doesn't include exogenous variables but I don't do it when it includes ones. > > If for example $p=q=1$ and $r=0$ the model is $$y_t = \alpha_0 + \alpha_1 > y_{t-1} + \tau_t + \varepsilon_t,$$ Let assume that the $y_t$ process is mean stationary, then $$E(y_t) = \mu,$$ for all t. Then $$E(y_t) = \alpha_0 + \alpha_1 E(y_{t-1}) +E(\tau_t) \rightarrow \mu=\alpha_0 +\alpha_1\mu +E(\tau_t)\rightarrow \mu = \frac{\alpha_0 +E(\tau_t)}{1-\alpha_1}. $$ So the process is mean stationary if and only $\alpha_1 \neq1$ and the process $\tau_t$ is mean stationary. Could you please help me with this. # Question 2: I would like to ask why do we want a time series to be stationary / invertible. What is the intuition behind this? What are the consequences of non-stationarity. Would you please recommend any reference to clarify this theme? Also if I estimated the unknown parameters of the model without imposing restrictions to them, this would influence my results? In what sense?


title: Why does R's lm() return different coefficient estimates than my textbook?

text: ## Background I'm trying to understand the _first_ example in a course on fitting models (so this may seem ludicrously simple). I've done the calculations by hand and they match the example, but when I repeat them in R, the model coefficients are off. I thought the difference may be due to the textbook using population variance ($\sigma^2$) whereas R may be using sample variance ($S^2$), but I can't see where these are used in the calculations. For example, if `lm()` uses `var()` somewhere, the help section on `var()` notes: > The denominator n - 1 is used which gives an unbiased estimator of the > (co)variance for i.i.d. observations. I have looked at the code for both `lm()` and `lm.fit()` and neither make use of `var()`, but `lm.fit()` passes that data to compiled C code (`z <\- .Call(C_Cdqrls, x, y, tol, FALSE)`) which I don't have access to. ## Question Can anyone explain why R is giving different results? Even if there is a difference in using sample vs population variance, why do the coefficient estimates differ? ## Data Fit a line to predict shoe size from grade in school.               # model data     mod.dat <- read.table(         text = 'grade shoe                     1    1                     2    5                     4    9'         , header = T);          # mean     mod.mu  <- mean(mod.dat$shoe);     # variability      mod.var <- sum((mod.dat$shoe - mod.mu)^2)          # model coefficients from textbook     mod.m  <- 8/3;     mod.b  <- -1;          # predicted values  ( 1.666667 4.333333 9.666667 )     mod.man.pred       <- mod.dat$grade * mod.m + mod.b;     # residuals         ( -0.6666667  0.6666667 -0.6666667 )     mod.man.resid      <- (mod.dat$shoe - mod.man.pred)     # residual variance ( 1.333333 )     mod.man.unexpl.var <- sum(mod.man.resid^2);     # r^2               ( 0.9583333 )     mod.man.expl.var   <- 1 - mod.man.unexpl.var / mod.var;          # but lm() gives different results:     summary(lm(shoe ~ grade, data = mod.dat))                    Call:     lm(formula = shoe ~ grade, data = mod.dat)          Residuals:           1       2       3      -0.5714  0.8571 -0.2857           Coefficients:                 Estimate Std. Error t value Pr(>|t|)     (Intercept)  -1.0000     1.3093  -0.764    0.585     grade         2.5714     0.4949   5.196    0.121          Residual standard error: 1.069 on 1 degrees of freedom     Multiple R-squared:  0.9643,    Adjusted R-squared:  0.9286      F-statistic:    27 on 1 and 1 DF,  p-value: 0.121      ## Edit As Ben Bolker has shown, it looks like teachers make mistakes sometimes. It seems that R calculations are correct. Moral of the story: don't believe something just because a teacher says it is true. Verify it for yourself!


title: Chow test multiple variables

text: great site! I have collected monthly deals data for how new ships are being finance from 2005-2012. The data consists of "deals" made by different shipping companies and the amount raised in each deal. The deals can split into the following categories: Bank Loans, Bonds, IPOs and FOs. Bank loans has been the premier financing option, accounting for 70-80% of total capital raised pre-2008. I am therefore looking to do a time series analysis to see if there has been a structural break (chow test) in the use of bank loans, and that other forms of financing such as Bonds, IPOs and FOs has increased post the financial crisis (2008). My question is, will I have to conduct a chow test for each variable to prove the structural break or can I somehow incorporate all variables into one analysis? Thanks in advance!


title: Interpreting a residuals vs fitted plot and extracting points

text: I'm doing a multivariate linear regression with R, and i find myself with the following residuals vs fitted plot: ![plot](http://i.stack.imgur.com/1azOo.png) As you can see there is a very regular line of points that seems to follow a precise pattern. **My questions are:**   1. How do I interpret such a behavior, and what can I do to fix it?   2. Is there a way to isolate/extract those points? I'd like to take a look at them individually in my data set to see if by examining them I notice some patterns in the data. **Additional info:** My model is:               v.lm = lm(sqrt(v.stima$Y)~., data=v.stima)      Y is a count-variable (non-negative integer). I'm using sqrt because without it the plot has the typical "funnel" shape that indicates a non-homoscedastic error.