title: How to combine WEKA classifiers

text: I need to utilize two different classifier to get best classification results. Since, it seems that they complement each other (not sure I am not expert btw). ROC characteristics are given below (testing scheme is 10-fold cross validation): **MetaCost [0 8; 1 0] Alternating Decision Tree (ADTree)**                              TP Rate   FP Rate   Precision   Recall  F-Measure   ROC Area  Class                      0.973     0.62       0.119     0.973     0.212      0.696    YES                      0.38      0.027      0.994     0.38      0.55       0.696    NO     Weighted Avg.    0.427     0.074      0.925     0.427     0.523      0.696      **MetaCost [0 8; 1 0] Logistic**                              TP Rate   FP Rate   Precision   Recall  F-Measure   ROC Area  Class                      0.604     0.161      0.245     0.604     0.348      0.771    YES                      0.839     0.396      0.961     0.839     0.896      0.769    NO     Weighted Avg.    0.821     0.377      0.904     0.821     0.853      0.769      I have tried `Voting` but could not get desired results. Therefore, it is time to seek for expert help. Could you please advise me a solution, if possible? Thanks in advance. Also, as a reminder I am not an expert. EDIT: Best I can get:               Vote combines the probability distributions of these base learners:         weka.classifiers.meta.AdaBoostM1 -P 100 -S 1 -I 10 -W weka.classifiers.meta.MetaCost -- -cost-matrix "[0.0 8.0; 1.0 0.0]" -I 10 -P 100 -S 1 -W weka.classifiers.trees.ADTree -- -B 10 -E -3         weka.classifiers.meta.AdaBoostM1 -P 100 -S 1 -I 10 -W weka.classifiers.meta.MetaCost -- -cost-matrix "[0.0 8.0; 1.0 0.0]" -I 10 -P 100 -S 1 -W weka.classifiers.functions.Logistic -- -R 1.0E-8 -M -1     using the 'Product of Probabilities' combination rule                          TP Rate   FP Rate   Precision   Recall  F-Measure   ROC Area  Class                      0.706     0.204      0.231     0.706     0.348      0.825    YES                      0.796     0.294      0.969     0.796     0.874      0.825    NO     Weighted Avg.    0.789     0.287      0.91      0.789     0.832      0.825


title: Is it possible to have a valid regression with debt to income as the D.V and income as one of the I.V's?

text: I have been assigned the following question "On the basis of the Bank of England's 2013 survey of the Financial position of British households, given in the data file [households], examine the statistical and econometric relationships **between household debt-to-income ratios and the respondents background characteristics, in particular: social grade; age-group; income; saving and housing tenure**. Consider the importance of your results with respect to the impact on household debt-to-income ratios of: (a) an increase in interest rates; and (b) a fall in the rate of price inflation. My question is is it valid to regress D/Y=f(Y,...) since Y appears on both sides? If not then what would be the way to deal with the problem. I thought perhaps loging the regression so that ln(D) can be regressed on ln(Y) so that an implicit relationship between D/Y and Y can be obtained after the regression. Any help will be much appreciated.


title: LIBSVM ROC PLOt

text: my question is "What is command used to generate ROC PLOT in libsvm as well as for specificity , sensitivity and precision what is the command in libsvm " i Am using ubuntu


title: Plotting a heatmap given a dendrogram and a distance matrix in R

text: I have dendrogram and a distance matrix. I wish to compute a heatmap -- without re-doing the distance matrix and clustering. Is there a function in R that permits this?


title: Absorption probability in 1D RW with asymmetric step sizes, $ x<0 $

text: What is the probability of absorption at $ 0 $, as a function of position $ x $, for a 1D random walk (on $ \mathbb{Z} $) with asymmetric step sizes? For example, suppose that you can take two steps to the right with probability $ p $, or one step to the left with probability $ q=1-p $. I calculated the probabilities for $ x>0 $ via difference equations*, but I'm stuck for the case $ x<0 $. Obviously it matters whether $ x $ is even or odd; in the case that $ p=1 $ the probabilities should bounce between $ 0 $ (if $ x $ odd) and $ 1 $ (if $ x $ even). *For $ x>0 $, I found $ P(x)=\left(\frac{-p+\sqrt{4p-3p^{2}}}{2p}\right)^{x} $.


title: Why do percentage-based forecast error measures assume a meaningful zero?

text: I've seen this comment made in various textbooks and papers. For example, the online textbook by Rob J Hyndman and George Athanasopoulos states at http://otexts.com/fpp/2/5/ that > Another problem with percentage errors that is often overlooked is that > they assume a meaningful zero. For example, a percentage error makes no > sense when measuring the accuracy of temperature forecasts on the > Fahrenheit or Celsius scales. Why is this the case?


title: Recommended reading for understanding when the bootstrap will fail?

text: It is known that the bootstrap can fail. I read in Section 6 of Bickel and Freedman (1981) that the bootstrap fails when you wan to use it to evaluate the MLE for estimating the parameter of a continuous uniform distribution. I read Secion 7.4 of the book by Efron and Tibshirani but I'm not able to find the reference they pointed to. Could someone point me to some more easily accessible things that I could refer to? Thanks!


title: Smoothing when standard errors are known/estimated

text: I have a set of estimates of function values, along with estimates of their standard errors. To somewhat simplify matters: for $x$ running from $1$ up to $70$ in (integer) steps (in fact a parameter in model selection), I have an estimate of the matching $y$-value (in fact an AUC) and an estimate of its SE. I cannot make assumptions on the form of the curve (it's not even certain that it's unimodal, for example), other than that it should be smooth on the scale of my $x$ values. I am looking for a way to translate these to a more smooth form, because the method I am using tends to create the occasional spike that I hope to circumvent. Now, most smoothing algorithms only use the $y$ values, which, in this case, seems to be a waste of information (i.e. not using the SE). In addition, I would also like to smooth out an error band from the error flags. Somewhat to my surprise, I have not found a technique that provides this. Have I not looked well enough? I would also like to perform this in R, so if anybody has a package that does the trick, I'd be very happy.


title: How to learn fitting models?

text: * I am new to R(learned the basics)   * In my current job I have to solve the following problem(the following is an example)  >   * Suppose we have 10, 000 users >   * For each user we have his expected monthly spending(in dollars)(with > some method that I don't know) for one year >   * What we have to do it fit a model where if for that year we give two > dates(Jan-1, March 16), the model should tell us the expected spending(in > dollars) for this duration > Question: >   * How to go about and fit a model >   * How to see graphically that model fits correctly with the data points I > have >   * How to predict and get answer > Please assist me   Thank you


title: How to select a radial basis function?

text: Currently I am investigating interpolation of 3D data with radial basis functions (RBF) and I am wondering that there are quite a few families of such (see table1 here). However, I cannot find any discussion for engineers about the intuitive choice of them. The are basically these aspects which differ across the RBF families:   * truncated power function with polynomials slide 21   * oscillatory RBF slide 36   * Euclid's hat functions slide 44   * Buhmanns RBF slide 49 Looks pretty much like a polynomial RBF from Wendland to me :( So why do we need so many of them? In my work I have to choose one of them, but I cannot justify the choice. :( So if someone asks me for that, I just tell that "I had a good dice"...


title: Multilevel model: variance significant, but covariance is not. What to conclude?

text: I estimate a two-level model, random slopes and random intercepts. The variance of each random term is significant based on Chi-square test, but the covariance (of slopes and intercepts) is not significant. Can I just safely remove the covariance term, thus simplifying the variance function? Are there any substantive conclusions following from that? Thanks in advance.


title: Logistics regression decision categorical variable makes huge difference in Kaggle-Titanic problem in R

text: I have this linear model:               fit = glm(Survived ~ Pclass + Sex + Age + SibSp + Parch + Fare,                data=passengers, family=binomial)     summary(fit)           Coefficients:              Estimate Std. Error z value Pr(>|z|)         (Intercept)  5.318162   0.571693   9.302  < 2e-16 ***     Pclass      -1.175648   0.145979  -8.054 8.04e-16 ***     Sexmale     -2.760823   0.199952 -13.807  < 2e-16 ***     Age         -0.043866   0.008220  -5.336 9.49e-08 ***     SibSp       -0.428252   0.106963  -4.004 6.24e-05 ***     Parch       -0.099051   0.118328  -0.837    0.403         Fare         0.002587   0.002362   1.095    0.274       So according to the Z score I can say that the decision variables `Pclass`, `Sexmale`, `Age`, `SibSp` deal a huge part in deciding whether the depending variable will be 0 or 1. Now I change the fit to the following (added Title):               fit = glm(Survived ~ Pclass + Sex + Age + SibSp + Parch + Fare + Title,                data=passengers, family=binomial)     summary(fit)          Coefficients:               Estimate Std. Error z value Pr(>|z|)         (Intercept)  3.458e+01  1.906e+03   0.018  0.98552         Pclass      -1.108e+00  1.556e-01  -7.123 1.06e-12 ***     Sexmale     -3.180e+01  1.906e+03  -0.017  0.98669         Age         -3.110e-02  9.574e-03  -3.249  0.00116 **      SibSp       -6.079e-01  1.250e-01  -4.862 1.16e-06 ***     Parch       -3.614e-01  1.346e-01  -2.685  0.00726 **      Fare         4.168e-03  2.573e-03   1.620  0.10531         TitleDr     -6.257e-01  1.691e+00  -0.370  0.71142         TitleLady   -1.605e+01  1.452e+03  -0.011  0.99118         TitleMaster  2.541e+00  1.552e+00   1.637  0.10163         TitleMiss   -2.982e+01  1.906e+03  -0.016  0.98752         TitleMlle   -1.642e+01  2.356e+03  -0.007  0.99444         TitleMr     -9.856e-01  1.439e+00  -0.685  0.49339         TitleMrs    -2.904e+01  1.906e+03  -0.015  0.98784         TitleMs     -1.498e+01  3.064e+03  -0.005  0.99610         TitleRev    -1.576e+01  9.616e+02  -0.016  0.98692         TitleSir    -4.003e-01  1.705e+00  -0.235  0.81438       So now my model is much "weaker" in a sense that I have more weak decision variables and the `Sexmale` is also "overshadowed" with the categorical `Title` variable. Could someone tell me the reason? Also, please help me deciding if the new model is really weaker than the previous one or am I missing some fundamental thing here?


title: Consistent ranked list for ROC interpolation

text: For classifiers with binary outputs, their performance is summarized by a true positive rate and false positive rate. To interpolate the performance between two classifiers $A$ and $B$ with their respective (true positive rate, false positive rate) being $(t_A, f_A)$, and $(t_B,f_B)$, we can construct a new classifier to achieve true positive rate of $t_A+(t_A-t_B)\frac{f_C-f_B}{f_A- f_B}$ at a false positive rate of $f_C$ by sampling. For any data point $x$, take the label $\hat{y}_C$ from the classifier $A$ (i.e., $\hat{y}_C \leftarrow \hat{y}_A$) with the probability $p_A = \frac{f_C-f_B}{f_A-f_B}$, and from $B$ (i.e., $\hat{y}_C \leftarrow \hat{y}_B$) with probability $1-p_A$. My question is, can we interpolate between at points so we achieve a ranked list? i.e., assuming $f_B<f_A$, as we go from $f_B$ to $f_A$, we keep adding to the list of points predicted to be Class 1?


title: Which to believe: Kolmogorov-Smirnov test or Q-Q plot?

text: I'm trying to determine if my dataset of continuous data follows a gamma distribution with parameters shape $=$ 1.7 and rate $=$ 0.000063. The problem is when I use R to create a Q-Q plot of my dataset $x$ against the theoretical distribution gamma (1.7, 0.000063), I get a plot that shows that the empirical data roughly agrees with the gamma distribution. The same thing happens with the ECDF plot. However when I run a Kolmogorov-Smirnov test, it gives me an unreasonably small $p$-value of $<1\%$. Which should I choose to believe? The graphical output or the result from KS- test? ![QQplot and ECDF plot](http://i.stack.imgur.com/4UctD.png)


title: Standard deviations and confidence intervals (weighted) running average

text: My question is related to this one. I am calculating averages, actually as many as I have samples because I calculate a running average, and for equal weighting I know how to calculate the $95\%$ CI, with $\mu\pm1.96*\sigma_{original}$. However there are two things that I dont know how to do. Can I get the $\sigma_{original}$ from the same method (convolution through Fourier transform)? How do I calculate the $95\%$ CI if the weighting is not equal (so if I convolute with a window that is not flat)?


title: How to report SPSS OLS output?

text: Consider the model below: ![enter image description here](http://i.stack.imgur.com/1cs1n.png) In many research papers, significance of statistical results are indicated by `*`, `**`, and `***` as a significant value at respectively the `1%`, `5%`, or `10%` level. Does a `Sig.` value of `.047` (forninstance, `LAW_LE`) mean a significane at the 4,7% level? (So that I may report `LAW_LE 0.313**`)? Second, I see many times the Std.Error is also reported in papers. How do I need to interpret the standard error in combination with the significance? For instance, `COM_BIK` has a standard error of 0.045 and a significance of 0.021; whereas `HOF_MAS` has a low standard error (0.002) but a significance of .416. What can be said about these two variables based on this information? Your help is greatly appreciated.


title: Analysis of discrete (count) time-series data

text: I have a 2X2 factorial experiment where I am interested in seeing the effect of two different nutrient solutions (N and W) on the appearance of root tips of two different plant species (A and B). I used 3 plants per species X treatment combination. I took daily counts of the number of root tips of each plant over the course of 9 days. The data look like this. ![enter image description here](http://i.stack.imgur.com/ZS3SQ.png) I have never worked with discrete data before. I was thinking of fitting an exponential growth function such as: Tips=T0e^rt , where T0 is the number of tips at time 0 (the observed initial number of tips), r is the relative rate of addition of new tips and t is time (in days, considered continuous?), and see if there were differences in the r parameter between my species X treatment combinations (or if there are only main effects…). Any suggestions? I know I will have to account for issues such as correlated observations and having (random?) subject-specific effects on the model parameters. I have no idea how to do this in the context of discrete data analysis, and would really appreciate any suggestions of potential models to use, pitfalls and references that would give me some background.


title: Anyone know of a simple dendrogram visualizer?

text: I've written a small hierarchical clustering algorithm (for better or for worse). I'd like a quick way of visualizing it, any tooling ideas?


title: Is it legitimate to use a conditional PDF derived using kernel density estimation for hypothesis testing?

text: Suppose I have some sample $X$ drawn from some unknown multivariate distribution $F(A,B)$, and I want to test the null hypothesis that a particular point $x$ was drawn from $F$. Would it be legitimate to fit the PDF for $F$ using kernel density estimation, then evaluate it at $x$ and take this as my p-value? Likewise, if I'm interested in the conditional probability $p(x~|~A=a)$, could I basically use the same approach, but just evaluate the PDF conditioned on $A = a$?


title: How to do Empirical Bayes Estimation of HLM parameteres in R?

text: I am modelling the impacts on students' achievemens using HLM (2 levels only), with the lme4 package in R. I would like now to estimate the Empirical Bayes Estimates and the Empirical Bayes Grand Means for each group of the 2nd level, using R. Please, would you be able to point me to any source or give me any tips on that? Much apreciated.


title: Produce a list of variable name in a for loop, then assign values to them

text: I wonder if there is a simple way to produce a list of variables using a for loop, and give its value.               for(i in 1:3)     {       noquote(paste("a",i,sep=""))=i     }      In the above code, I try to create `a1`, `a2`, `a3`, which assign to the values of 1, 2, 3. However, R gives an error message. Thanks for your help.


title: AIC values and their use in stepwise model selection for a simple linear regression

text: The Wikipedia article for AIC says the following (emphasis added): > As an example, suppose that there were three models in the candidate set, > **with AIC values 100, 102, and 110.** Then the second model is > exp((100−102)/2) = 0.368 times as probable as the first model to minimize > the information loss, and the third model is exp((100−110)/2) = 0.007 times > as probable as the first model to minimize the information loss. > > **In this example, we would omit the third model from further > consideration.** We then have three options: (1) we could decide to gather > more data, in the hope that this will allow clearly distinguishing between > the first two models; (2) we could simply conclude that the data is > insufficient to support selecting one model from among the first two; (3) we > could take a weighted average of the first two models, with weights 1 and > 0.368, respectively, and then do statistical inference based on the weighted > multimodel. However, a video discussing the stepwise method for model selection in R removes the **smallest** AIC value . It may be that I am grossly misunderstanding something in between how AIC works and how AIC is applied. Could anyone explain why we would not want to select the largest value in the video as was done in the Wikipedia example?


title: Variable selection and forecasting with regARIMA models

text: I have a couple of questions about regARIMA models:   1. What is the underlying principle of the R function `auto.arima` when xreg is different from NULL? Does it first perform a regression of the time series on the explanatory variables and then select the best ARIMA model for the residuals of the previous regression? If I am not wrong, this would mean that the coefficients of the regression would be the same whether or not we say that the error terms are i.i.d normal variables with 0 mean and constant variance or follow an ARIMA process (in other words, whether or not you fit a simple regression model or a regARIMA model). Thus, fitting a regARIMA model would not have any impact on the coefficient of the explanatory variables...   2. I read somewhere (maybe on this forum) that it was not relevant to look at the t-statistics or p-values to determine which variables have to be included in the model (in particular the dummy variables representing the potential seasonalities). I do not manage to find the link again. Could anyone remind me why?   3. I did not manage to find any details about forecasting regARIMA models. I would say that the predicted value of your time series at (say) t+h is the sum of the fitted value given by the regression model you have fitted (at t+h) plus the predicted value of your error term at t+h (obtained via the usual forecasting methods for ARIMA models). Thank you in advance for your precious help.


title: Selecting regression model for a non-negative integer response

text: I have a series of non-negative integers $y=(y_1,y_2,..., y_n)$ and a design matrix $y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1 x_2$, where $x_0$ and $x_1$ are $0$ or $1$, $x_1x_2$ is the interaction, and $\beta_0 \ldots \beta_3$ are parameters we want to estimate. For example, the data look like               y    x1    x2    x1*x2     10   0     0     0     23   0     1     0     18   1     1     1     19   1     0     0     25   0     1     0     ...      I want to estimate the $\beta_0$, $\beta_1$, $\beta_2$ and $\beta_3$ coefficients and perform a test to see if any coefficient is nonzero. **There are several different regression models** that might be applied to this case:   1. Simple linear regression: `lm`   2. Poisson regression (when $y$ follows a Poisson distribution): `glm` with family = poisson   3. Quasi-poisson regression (when $y$ is over-dispersed; that means $\text{sd}(y) \gt \text{mean}(y)$): `glm` with family = quasi-poisson   4. Negative binomial regression (when $y$ is over-dispersed, $\text{sd}(y) \gt \text{mean}(y)$): `glm.nb`, in MASS package. **The questions I want to ask are:**   1. How should I select the model for this dataset? Is there any way to choose the right model based on some descriptive statistics of my dataset?   2. How should I check and validate if the fitted selected model is right for my data? Thanks everyone.


title: Tree structured proportion of success : statistical significance?

text: I am exploring a tree structured as follows: A/ parent node = parent population, size N, with a proportion of success for an experiment (eg. p = nb success / N) B/ child nodes = sub populations for a give parent node. Ex.: level one child sizes are n1, n2, ... nn, where n1 + n2 + ... + nn = N. Success proportions are also computed (p1 = nb1 success / n1, ..., pn = nbn success / nn) This could be summarized by the enclosed picture. ![The problem](http://i.stack.imgur.com/6DJXT.png) What I am trying to do is cut the nodes where the proportions of success are not statistically significant. In other words, how can I trust the proportion of success for a given node, according to its population size? The context of this question is pretty specific: p, p1, ... are always very low (between 0 and 0.0002, and node size can vary a lot, from 10 to several thousands of observations. So how can I compute the minimum node size, so that I can consider a node proportion of success as statistically significant? Moreover, do I have to do a computation for each of the tree node, do I have to take the dependencies between the nodes into account, or is a global approximation for all the nodes acceptable? I'm working with R. My first insight was use the prop.test function for all node and check the p-value, but I realize this completely wrong. As anyone an idea to deal with this tricky problem?


title: How to solve Chuck-a-Luck puzzle

text: ### Story (spoilers!) The puzzle: > Chuck-a-Luck is a gambling game often played at carnivals and gambling > houses. A player may bet on any one of the numbers 1,2,3,4,5,6. Three dice > are rolled. If the player's number appears on one, two, or three of the > dice, he receives respectively one, two, or three times his original stake > plus his own money back; otherwise he loses his stake. What is the player's > expected loss per unit stake? (Actually the player may distribute stakes on > several numbers, but each such stake can be regarded as a separate bet.) You can find a solution in Fifty challenging problems in probability with solutions (F. Mosteller). Before reading the solution of course, I worked on my own, and I thought "hey, there is nice symmetry here, so let's assume I bet on 1 and then calculate loss". The Author thought (just an assumption) "hey, there is nice symmetry here, so let's assume I bet at the same time on all numbers". I calculated loss about 11%, the Author calculated loss at ~8%. I rechecked my calculations, they seem fine, I rechecked the Author's calculations and they also look fine. However, no matter how long I stare at those outcomes, 11% is not 8%. _Btw. you can find the same comparison one by one here:http://www2.washjeff.edu/users/mwoltermann/Mosteller/contents.htm -- number 6._ This made me wonder... ### Question Are those approaches actually both incorrect? I mean, the question is plain and simple **"What is the player's expected loss per unit stake?"**. You can model the scenarios, but adding something to the question like "... assuming the player bet on all the numbers" **changes** the question. So Author answers on similar, yet different question. So, **MY** question is -- am I right with above suspicion, that in order to correctly solve this problem one should calculate expected loss per unit stake period. The calculation should be done really not with parameter (bet on 1) or (bet on all) but (bet A on 1, bet B on 2, bet ..., bet F on 6) and only with such parameters calculate the expected loss. Am I right with my new approach?


title: How do you find the average in this Poisson Distribution problem?

text: Suppose that the number of accidents occurring at a particular intersection each week is a Poisson random variable. It is estimated that the probability of having at least one accident in a week, at the intersection, is 0.98. Calculate the average number of accidents occurring each week at the intersection. Thanks in advance for the help!


title: Simulation of an exponentially distributed variable

text: A one-period jump diffusion model of a stock is the following: given the stock value at $T0$, $S(T0)$, the stock value at time $T$, $S(T) = S(T0)exp(x)exp(y)$. Given the inputs $m$, $s$, $pup$, $pdown$ and $\lambda$, and:   1. x is normally distributed with mean m and standard deviation s   2. y equals 0 with probability 1 – pup – pdown (i.e. no jump).    3. y is positive with probability pup; if it is positive , y is exponentially distributed with density p(y) = λexp(-λy) for y > 0.   4. y is negative with probability pdown; if it is negative, y is exponentially distributed with p(y) = λexp(λy) for y < 0.   5. x and y are uncorrelated & independent. I am to write a MATLAB program that siumulates the stock prices. Two random numbers per trial are required:   1. A uniformly distributed random number that determines if a jump has occurred   2. An exponentially distributed random number that determines the size of the jump I've transormed a uniform random variable to an exponential variable using $-ln(1 - U(a, b)) / \lambda$. Where I'm confused is how we're using the uniformly distributed random number that determines if a jump has occurred given I am given pup and pdown. Also, why do I need an exponentially distributed random number for the jump size. I thought I just integrate the pdf? For some reason I struggle mightily with the exponential distribution.


title: libsvm grid search on Parallel Computing Toolbox Matlab

text: I am using the following code for grid search on libsvm:               bestcv = 0;     for log2c = -1:3,     for log2g = -4:1,        cmd = ['-v 5 -c ', num2str(2^log2c), ' -g ', num2str(2^log2g)];        cv = svmtrain(heart_scale_label, heart_scale_inst, cmd);        if (cv >= bestcv),          bestcv = cv; bestc = 2^log2c; bestg = 2^log2g;        end       fprintf('%g %g %g (best c=%g, g=%g, rate=%g)\n', log2c, log2g, cv, bestc, bestg, bestcv);      end     end      To increase the speed, I can either go for Parallel computing toolbox or running this on GPUs (though I guess first option would be better because GPUs are good at operations like convolution etc which are not so heavy in this case - please correct me if I am wrong.) How can I optimize this by running it on Parallel computing toolbox? I am beginner in this field and therefore need some input on how to implement the parfor loops. **EDIT** Can it be done using matlabpool?


title: Testing if code is generating random variables consistent with theory

text: In software development, we often write unit tests to ensure that a particular portion of the code is working as expected. For example, consider the problem of drawing a color at random from Colors=(Red, Green, Blue). Suppose, that we write a function:               drawn_color = draw_color(Colors)      One way to assess whether the code is working correctly or not is to do the following:               tot_iter = 10000       For iter = 1 to tot_iter          drawn_color(iter) = draw_color(Colors)       Next       And then check if the frequencies of the colors in drawn_color are roughly 10000/3. However, the above assessment requires visual inspection and I would rather automate the whole thing so that I can write a test code like so:               If abs(empirical freq - expected freq) < = epsilon then the code is valid.      Two questions:   1. How does epsilon vary as a function of tot_iter?   2. Is there a generic approach to above the issue where the random variables could be drawn from other distributions/densities (e.g., normal, uniform, discrete etc)?


title: New Systematic Quantile Normalization

text: I think I've developed a systematic quantile normalization technique. I did this with music but I think it can also be done with light and other frequency based information. The algorithm is as follows: Given a series of frequency data: >   1. Systematically move it to a known subspace containing only periodic and > rounded samples of the data at that frequency. That is cause the data to > inhabit the same vector space but not yet in the same Cartesian space. This > means transcribe the frequencies into a system of music or simular > representation. Western music is a periodic sample expressing a frequency > but rounded to a few Hz. > >   2. Transpose the frequency data into the same range of frequency. Moving > average range from 300 Hz to 1000 Hz to match the range of the one frequency > of another sample. This is move to a similar Cartesian space. This would be > moving Soprano, Alto, Tenor and Baritone voices into alignment with each > other. > >   3. Normalize the frequency data into the exact same range of frequency > inducing loss of data. Points which unbalance the range are moved up higher > or lower. Attempt to overlap in Cartesian space. This is making the vocal > data have the same approximate same "root" and "key". > >   4. Correlate local maxes and min at similar percentile spacing to > establish levels of similarity. After transcribing the music, moving the > range and key: ignore the spurious notes generated in between major notes > from background noise. Since no two musical performance are exactly time > synchronized a percentage can be used to identify when notes are at a > simular place. > > 0 Frequency ![freq](http://i.stack.imgur.com/dXIFt.png) 1 Move to system ![movetomusicsystem](http://i.stack.imgur.com/cEY82.png) 2 Transpose normalize ![transpose](http://i.stack.imgur.com/q6nQW.png) 3 Correlate ![corrolate](http://i.stack.imgur.com/IdME9.png) This is a systematic quantile extrapolation into a common vector space denoted by range and tone based on a loss of data when points are at extreme or highly variable beyond the ability of the subspaces rounding and sampling. I have proven to my satisfaction that it is able to recognize a pattern of similarity in tone between any two frequencies at different pitches. Could someone verify or disprove (and better qualify) this algorithm as quantile (in normalization type) and differential (in pitch tone comparative determination) please?


title: Can I trust a model if I do not check the prediction error?

text: My main field is machine learning and 90 % of what I do is to try to improve the prediction error. Recently I have started to work with a medical group. They are mainly doctors so I do not know how much I can trust their statistical knowledge. What they commonly do is to fit logistic regression and without any consideration of the model performance they start to look to the coefficient in order to understand which one is more important or stuff like that. In a similar way I fitted a random forest model and the I have done the importance plot something like this even if the predicted performance was very poor. How much can I trust these kind of results in model where the prediction performance is poor?


title: Observer agreement with Cohen's Kappa - is there a way to combine ratings to investigate if one rating impacts the other

text: I am a medical scientist doing observer agreement analysis on some very simple data using Cohen's Kappa. I have 2 observers. They score a bunch of scans on a 3 point scale for a pattern called honeycombing (definitely present, possibly present, absent). They also score for emphysema on the same 3 point scale. It is commonly believed that observers confuse these patterns. The unweighted kappa for honeycombing is 0.56. The unweighted kappa for emphysema is 0.17. I need to demonstrate if the emphysema pattern does confuse observers when trying to identify honeycombing and I assumed that if I combine the scores in some way, if this hypothesis is true, the kappa for honeycombing should improve. - i.e combine the scores so that the result represents the observers determination for the present of honeycombing OR emphysema. I am confused however - when I create a binary score for honeycombing (present/possibly present versus absent) and the same for emphysema and then combine these binary scores to give a "honeycombing or emphysema present or absent" the kappa falls to 0.46. I find this hard to explain. I am a novice to statistics. Opinions appreciated - I am sure its a point of logic I am missing.


title: Repeated Measures ANOVA with several measures entered together: how to find main effect

text: I've run a 2 (treatment and no treatment) x 2 testing occasions (pre and post) Repeated Measures ANOVA. In addition, I entered several tests given at both time points, so under the Repeated Measures Define Factor(s) window option for Measure Name, I've entered 5 different tests used. In the output, I am trying to determine which table indicates the main effect _across_ all tests. What I am seeing instead is a Tests of Within-Subjects Effects Multivariate table (that I am ignoring) and a Univariate table that lists `session`, `session*group`, and `error(session)` effects by test. To report a main effect to address the question of whether _any_ tests differed between the groups from pre to post, where would I find that information?


title: Does stepwise regression provide a biased estimate of population r-square?

text: In psychology and other fields a form of stepwise regression is often employed that involves the following:   1. Look at remaining predictors (there are none in the model at first) and identify the predictor that results in the largest r-square change;    2. If the p-value of the r-square change is less than alpha (typically .05), then include that predictor and go back to step 1, otherwise stop. For example, see this procedure in SPSS. The procedure is routinely critiqued for a wide range of reasons (see this discussion on the Stata website with references). In particular, the Stata website summarises several comments by Frank Harrell. I'm interested in the claim: > [stepwise regression] yields R-squared values that are badly biased to be > high. Specifically, some of my current research focuses on estimating population r-square. By population r-square I refer to the percentage of variance explained by the population data generating equation in the population. Much of the existing literature I am reviewing has used stepwise regression procedures and I want to know whether the estimates provided are biased and if so by how much. In particular, a typical study would have 30 predictors, n = 200, alpha of entry of .05, and r-square estimates around .50. What I do know:   * Asymptotically, any predictor with a non-zero coefficient would be a statistically significant predictor, and r-square would equal adjusted r-square. Thus, asymptotically stepwise regression should estimate the true regression equation and the true population r-square.   * With smaller sample sizes, the possible omission of some predictors will result in a smaller r-square than had all predictors been included in the model. But also the usual bias of r-square to sample data would increase the r-square. Thus, my naive thought is that potentially, these two opposing forces could under certain conditions result in an unbiased r-square. And more generally, the direction of the bias would be contingent on various features of the data and the alpha inclusion criteria.   * Setting a more stringent alpha inclusion criterion (e.g., .01, .001, etc.) should lower expected estimated r-square because the probability of including any predictor in any generation of the data will be less.   * In general, r-square is an upwardly biased estimate of population r-square and the degree of this bias increases with more predictors and smaller sample sizes. ### Question So finally, my question:   * **To what extent does the r-square from stepwise regression result in a biased estimate of population r-square?**   * **To what extent is this bias related to sample size, number of predictors, alpha inclusion criterion or properties of the data?**   * **Are there any references on this topic?**


title: How to analyze paired and longitudinal data

text: I want your advice regarding data I have, that comes from the Ophthalmic world. I have a small sample of 12 people (in future studies will be larger, this is a small pilot). Each person received a medicine to lower his intraocular pressure, which is a metric scale measure (continuous). In this small study, prior giving the medication, the pressure was measured in both eyes to obtain baseline recordings. After giving the medication, the pressure was measured again in both eyes, at 5 additional follow-up visits of the patients in the clinic. Prior the intervention, the difference between the eyes was 0 on average among all 12 patients. The reason why the controlled eye was used in the first place, is because there is a large variability of the pressure with time, without any relation to the treatment, and there is a correlation between the variability in both eyes, only where the medication was given, it is assumed the entire line will be lower. I want to know for how long the medication takes effect. If I work with the differences, which is what my logic is, I am in big trouble, since after baseline, there is a fall (the difference is negative, I take D = treatment - control), then it start to go up, but not always linearly. I can't do a nonlinear mixed model of paired differences with 12 observations... I wanted to ask you how would you analyze such data ? I know how to analyze paired data (paired t-test), and I know how to analyze longitudinal data (GEE), however, how do you combine both in one model ? Thank you !


title: Covariance of ARMA(2,1) series

text: Consider the ARMA(2,1) time series $$ x_t−0.1x_{t−1}−0.06x_{t−2}=w_t−0.5w_{t−1} , $$ where $w_t$ is white noise with mean zero and variance $\sigma^2_w$. Find the expectation of $$ (x_t−0.1x_{t−1}−0.06x_{t−2})x_{t−h}, \quad h=0,1,... $$ and $$ (w_t−0.5w_{t−1})x_{t−h}, \quad h=0,1,... $$ and show how this can be used to formulate a difference equation to determine the covariance $γ(h)=E[x_{t+h}x_t]$. I'm using the book 'Time Series Analysis and Its Applications' Third Edition by Shumway and Stoffer. But can't find the relevant examples to draw on.


title: Standard error of slopes in piecewise linear regression with known breakpoints

text: **The situation** I have a dataset with one dependent $y$ and one independent variable $x$. I want to fit a continuous piecewise linear regression with $k$ known/fixed breakpoints occurring at $(a_{1}, a_{2}, \ldots, a_{k})$. The breakpoins are known without uncertainty, so I don't want to estimate them. Then I fit a regression (OLS) of the form $$ y_{i} = \beta_{0} + \beta_{1}x_{i} + \beta_{2}\operatorname{max}(x_{i}-a_{1},0) + \beta_{3}\operatorname{max}(x_{i}-a_{2},0) +\ldots+ \beta_{k+1}\operatorname{max}(x_{i}-a_{k},0) +\epsilon_{i} $$ Here is an example in `R`               set.seed(123)     x <- c(1:10, 13:22)     y <- numeric(20)     y[1:10] <- 20:11 + rnorm(10, 0, 1.5)     y[11:20] <- seq(11, 15, len=10) + rnorm(10, 0, 2)      Let's assume that the breakpoint $k_1$ occurs at $9.6$:               mod <- lm(y~x+I(pmax(x-9.6, 0)))     summary(mod)          Coefficients:                         Estimate Std. Error t value Pr(>|t|)         (Intercept)          21.7057     1.1726  18.511 1.06e-12 ***     x                    -1.1003     0.1788  -6.155 1.06e-05 ***     I(pmax(x - 9.6, 0))   1.3760     0.2688   5.120 8.54e-05 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      The intercept and slope of the two segments are: $21.7$ and $-1.1$ for the first and $8.5$ and $0.27$ for the second, respectively. ![Breakpoint](http://i.stack.imgur.com/2GsoE.png) * * * **Questions**   1. How to easily calculate the intercept and slope of each segment? Can the model be reparemetrized to do this in one calculation?   2. How to calculate the standard error of each slope of each segment?   3. How to test whether two adjacent slopes have the same slopes (i.e. whether the breakpoint can be omitted)?


title: A method to find the inconsistency or variation in the data

text: I am running an experiment (it's an image processing experiment) in which I have a set of paper samples, and each sample has a set of lines. For each line in the paper sample, its strength is calculated, which is denoted by say ' **s** '. For a given paper sample I have to find the variation amongst the strength values ' **s** '. If the variation is above a certain limit, we have to discard that paper. 1) I started with the standard deviation of the values, but the problem I am facing is that for each sample, order of magnitude for **s** (because of various properties of line like its length, sharpness, darkness etc.) might differ and also the calculated standard deviations values are also differing a lot in magnitude. So I can't really use this method for different samples. Is there any way where I can find that suitable **limit** which can be applicable for all samples. I am thinking that since I don't have any history of how the strength value should behave (for a given sample depending on the order of magnitude of the strength value more variation could be tolerated in that sample whereas because the magnitude is less in another sample, there should be less variation in that sample), I first need to find a way of baselining the variation in different samples. I don't know what approaches I could try to get started. Please note that I have to tell variation between lines within a sample whereas the **limit** should be applicable for any good sample.


title: difference in pre and post populations

text: I'm looking at the difference in demographics of people with a disease pre- outbreak and post-outbreak (e.g. age, sex, location etc). I know how to use $t$ tests and $\chi^{2}$-tests for comparison (e.g. difference of means for age), and I'm getting ridiculous statistical significance, but what I would like to know is the difference at the individual group level, i.e. at a state level, was there a difference pre/post outbreak? I have a column with "state" and a column with "outbreak", so how could I perform that analysis? I'm using Stata 13 (as per the wishes of my supervisors). If this question has already been answered could you please direct me there? **Edit** I have the seven Australian states (e.g. WA, NSW etc.) and I would like to know if there was a difference in the number of notifications pre/post outbreak, or if there was a difference in age pre/post outbreak for each state. Does that help explain the "individual group level"?


title: Garch model detection

text: Does the plotted data look like a garch model? If it does, with which test, I can support this claim ![http://i.stack.imgur.com/gEHr3.png](http://i.stack.imgur.com/gEHr3.png)               [1]   0.187415738   0.477868164   0.295777631   0.575167970   0.469157108       [6]   0.328171826   1.292894322  -0.847711535   0.392077876   0.116540719      [11]   0.314835442  -0.131731938   0.090753214   0.326724648  -0.504309863      [16]   0.522419894  -0.458490815  -0.922188177   1.312544072  -0.124784861      [21]  -0.440201131  -0.376074075   0.484759007  -0.544084907   0.184140601      [26]   0.634605831  -1.015191304  -0.392407326   0.069140836   0.837570607      [31]  -0.734046990   0.009650584  -0.297811017  -0.894844413   1.254269251      [36]  -2.202675529   3.181080668   2.738269460  -1.531934991  -0.441972078      [41]  -0.111205739   0.347052040   0.812613401   1.056684583  -0.498294836      [46]  -0.318961547  -0.109851318  -0.118844912   0.404173405  -0.910582285      [51]  -0.353860467  -0.750922999  -0.703591979  -0.637060347  -0.282754755      [56]  -0.953784737   1.576556620   4.151024359  -5.425327217  -6.103292196      [61]   7.943888186 -13.463531790   5.869748465 -14.205872960  -5.045698842      [66]   6.001954472   4.930371803  -1.892510511   1.850741003  -0.954016972      [71]  -0.077456198   0.040999835  -0.446823649   0.127273533   0.127966102      [76]   0.216028033   0.626917525   0.282341328   0.328576480  -0.157644004      [81]  -0.937382208   1.287429037   0.900623061  -0.677630409   0.427116916      [86]   0.246673589   0.204713960  -0.177800003   0.290096733   0.758407267      [91]  -0.328556185  -0.330455346  -0.526235619  -0.489013356   0.475158246      [96]   0.430567333  -0.775044149   0.139062847   0.374927797  -0.528934679     [101]   0.217091044  -0.019686328   0.690038054  -0.450631659   0.530362959     [106]  -0.466588401   0.895308589   0.013337137  -0.257514294  -0.346425256     [111]   0.002790152   0.565433354  -1.038591166   1.439547986  -0.475994171     [116]  -0.573993908   0.844415185  -0.778530803  -0.743884821   1.258147231     [121]  -0.088087152   0.508910627  -0.071854449   0.410070978   0.471372207     [126]   0.094210278   1.334666193  -0.010790060  -0.305813138  -0.347626733     [131]  -0.467378439   0.372987762   0.481824018   0.415398181  -0.180618112     [136]   1.013615048  -0.156073342   0.356076568   0.169108466   0.275734683     [141]  -0.648042994   0.269037495   0.669168433  -0.239254253   0.796882449     [146]   1.336804388  -0.354919126  -0.525481668  -0.066679296  -0.645930518     [151]   0.415104329   1.070218821   2.305348685  -3.405189139  -3.799563891     [156]   0.721310801  -0.864861449  -1.169885359   0.052965774  -0.681096223     [161]  -0.227592654   0.538791441   0.123075072   0.471443863  -0.598493272     [166]  -0.575874729  -0.497991275  -1.149177081  -0.391882551  -0.724743988     [171]   0.151001413  -0.596898412   0.762378302  -1.740324011   0.255643898     [176]   0.417552215  -0.641130460  -0.775619616   0.968210920  -0.110971587     [181]  -0.143987039  -0.041898461   0.384757277   0.303919564   0.367943147     [186]  -0.123326308   0.506591242   0.631918157  -0.265859230   0.325193212     [191]   0.532946942   0.633710920   0.590115624   0.307122036   0.391484459     [196]  -0.036289366  -0.614109628  -0.638930907   0.881656376  -1.584292446                         Call:     garch(x = data2, order = c(0, 1))          Model:     GARCH(0,1)          Residuals:          Min       1Q   Median       3Q      Max      -2.57713 -0.68644  0.01396  0.79931  2.85119           Coefficient(s):         Estimate  Std. Error  t value Pr(>|t|)         a0   0.21135     0.04085    5.174 2.30e-07 ***     a1   0.87577     0.16163    5.418 6.02e-08 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Diagnostic Tests:         Jarque Bera Test          data:  Residuals     X-squared = 1.4612, df = 2, p-value = 0.4816                   Box-Ljung test          data:  Squared.Residuals     X-squared = 0.0027, df = 1, p-value = 0.9584               Call:     garch(x = data2, order = c(1, 1))          Model:     GARCH(1,1)          Residuals:           Min        1Q    Median        3Q       Max      -2.257745 -0.252623  0.006404  0.248729  1.680547           Coefficient(s):         Estimate  Std. Error  t value Pr(>|t|)         a0 3.497e+00   8.228e-01    4.251 2.13e-05 ***     a1 1.048e+00   8.054e-01    1.301    0.193         b1 2.371e-14   1.779e-01    0.000    1.000         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Diagnostic Tests:         Jarque Bera Test          data:  Residuals     X-squared = 143.2293, df = 2, p-value < 2.2e-16                   Box-Ljung test          data:  Squared.Residuals     X-squared = 15.6938, df = 1, p-value = 7.447e-05          Call:     garch(x = data2, order = c(1, 2))          Model:     GARCH(1,2)          Residuals:           Min        1Q    Median        3Q       Max      -3.301812 -0.268550  0.006167  0.257807  2.279465           Coefficient(s):         Estimate  Std. Error  t value Pr(>|t|)       a0 3.327e+00   2.675e+00    1.244    0.214       a1 1.341e-01   7.183e-02    1.866    0.062 .     a2 1.299e-01   1.559e-01    0.833    0.405       b1 1.113e-13   7.818e-01    0.000    1.000       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Diagnostic Tests:         Jarque Bera Test          data:  Residuals     X-squared = 470.4355, df = 2, p-value < 2.2e-16                   Box-Ljung test          data:  Squared.Residuals     X-squared = 60.2426, df = 1, p-value = 8.438e-15


title: (Spatial) distance between cluster means

text: I'd like to cluster points based on a distance criteria. As I want to cluster spatial points I am using euclidean distance and a hierachical cluster approach. In a final step I'd like to cut the dendrogramm at a specific distance to make sure that all cluster means are spatially separated by a minimum distance specified. Here the R code so far:               # Creating the spatial points     set.seed(1)     x <- runif(100,0,150)     set.seed(5)     y <- runif(100,0,150)     df <- data.frame(x,y)     plot(df)          # Clustering using a hierachical approach based on euclidean     hc <- hclust(dist(df,method = "euclidean"), method="average")          # Specify a distance of minimum 60 units between cluster means?     df$memb <- cutree(hc, h = 60)          # plot clusters in colors     plot(df$x,df$y,col=df$memb)          # calculate cluster means     cent <- NULL     for(k in 1:length(unique(df$memb))){           cent <- rbind(cent, colMeans(df[df$memb == k, , drop = FALSE]))     }     cent <- as.data.frame(cent)          # plot points on top of clusters     points(cent$x,cent$y,pch=15)      However when I am cross-checking the distance between the cluster means, I observe minimum distances that are slightly below my threshold (e.g. 60):               # Summary of distance between cluster means     summary(dist(cent))      I guess that is related to the clustering method "average" which is different from the cluster means I am calculating thereafter. However, how can I achieve, that the distance between cluster means (or centroids, etc) is at least 60 units? What should I use as a measure of cluster mean, and what can I use as a clustering method?


title: Posterior distribution of a random variable

text: Im not understanding the following; suppose $y \sim N (\mu,\sigma^2)$ and we have a prior $\mu \sim N (\mu_0, \sigma^2_1)$ Then we can figure out the posterior distribution. What i dont understand is suppose $y=\mu+e$ where $e \sim N(0,\sigma^2_e)$ how do I know that $y \sim N (\mu,\sigma^2)$ i.e. how do I know that Y is centered around $\mu$? Won't it be centered around $\mu_0$? Because; $E(y)=E(\mu)+E[e]=\mu_0$


title: Correlation between independent variables in Generalized linear model

text: short intro: I wish to explain seedling abundance in a savanna ecosystem with two factors. one factor is 'canopy' which is either 0 (under a tree) or 1 (not under a tree). the other factor is distance from forest edge for which i have 5 different distances. I also have three covariates :'tree density' (which was measured for every distance), grazingpressure (1 per site), herbecaous biomass (measured for every plot). I used a generalized linear model because I have count data and I have a lot of zero outcomes. The problem is that some of these independents are correlated with eachother: Distance and herb biomass (0,288), Distance and tree density (0,492) and Grazingpressure and herbiomass (0,360). I know I cannot keep them in my model like this, but I also dont want to exclude some of them. Grazingpressure should be removed ofcourse because the causation is clear. But i do feel that 'distance' from forest edge and 'tree density' explain 'seedling abundance' in a ecological different way and thats why I want to keep both. question: How do i make a accurate estimation of the main effects of the independent variables in a generalized linear model when they are correlated?


title: Friedman/Wilcoxon test. How to enter data for multiple repetitions?

text: Id like to run a non parametric test such as the Friedman test on a dependent variable of my 8 subjects who have done a task under four different conditions. I only found examples of how to input the data which had one value per condition for each subject: four vectors of 8 values (one for each subject) for each of the four conditions. 1 X X X X X X X X 2 X X X X X X X X 3 X X X X X X X X 4 X X X X X X X X However my 8 subjects have done the task not once but 16 times under each of these four conditions. How would I enter the data into the test such as in SPSS? Can I just append the 16 values for each subject to these lists, creating four vectors of 8x16=128 values? Or would I have to average the 16 values for each subjects beforehand to yield four vectors with 8 values each? Thank you!


title: Most significantly frequent category

text: I am measuring the frequency of consonants in a sample lexicon of a language. The results are as follows:               +-----+-----+-----+-----+----+----+-------+     |  p  |  t  |  k  |  m  | n  | r  | Total |     +-----+-----+-----+-----+----+----+-------+     | 278 | 256 | 122 | 189 | 83 | 72 |  1000 |     +-----+-----+-----+-----+----+----+-------+      It is easy to see that the frequencies are rather different, and $\chi^2$ test for goodness of fit confirms this ($\chi^2$ = 208.5345, df = 5, p < 0.001). The question I am actually really interested is the following: The most frequent consonant is _p_ , but is this significant? In other words, can I claim that _p_ is the most frequent consonant? Is _p_ significantly more frequent than the second-most frequent consonant _t_ , or could the difference between the two due to chance? What test should I use? Is it ok to simply use a $\chi^2$ test for goodness of fit on the frequencies of _p_ and _t_ only?


title: Comparing models using anova() function in R

text: From the documentation for `anova()`: > When given a sequence of objects, ‘anova’ tests the models against one > another in the order specified... What does it mean to test the models against one another? And why does the order matter? Here is an example from the GenABEL tutorial:                   >  modelAdd = lm(qt~as.numeric(snp1))         >  modelDom = lm(qt~I(as.numeric(snp1)>=2))         >  modelRec = lm(qt~I(as.numeric(snp1)>=3))          anova(modelAdd, modelGen, test="Chisq")         Analysis of Variance Table              Model 1: qt ~ as.numeric(snp1)         Model 2: qt ~ snp1           Res.Df  RSS Df Sum of Sq Pr(>Chi)         1   2372 2320                               2   2371 2320  1    0.0489     0.82          anova(modelDom, modelGen, test="Chisq")         Analysis of Variance Table              Model 1: qt ~ I(as.numeric(snp1) >= 2)         Model 2: qt ~ snp1           Res.Df  RSS Df Sum of Sq Pr(>Chi)         1   2372 2322                               2   2371 2320  1      1.77     0.18          anova(modelRec, modelGen, test="Chisq")         Analysis of Variance Table              Model 1: qt ~ I(as.numeric(snp1) >= 3)         Model 2: qt ~ snp1           Res.Df  RSS Df Sum of Sq Pr(>Chi)           1   2372 2324                                 2   2371 2320  1      3.53    0.057 .         ---         Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       How do I interpret these output? Thanks!


title: Are all sequences of of random (uniform) numbers also uniformly distributed?

text: If I take some sequences of random numbers generated by a random number generator with uniform distribution, will the resulting sequences be uniformly distributed as well? By example, if I have a generator that returns 1, 2, or 3, what are the probabilities to get [1, 1, 1] and [1, 2, 3]?


title: Estimation of competitiveness of firm

text: Recently I have taken up one assignment wherein I am mesuring the competitiveness of indian apparel export firm. To compute competitiveness I used a survey where firm owners were given 12 questions to rate their competency from 1 to 10. After response collection, I did a Factor Analysis (PCA) which gave me three factors. Besides I have also asked exporters to rate the macro economic factors and challenges they are facing in exports on a scale ranging from 1 to 5. Moreover structural variables like location, type of firm, employee strength, export performance, ownership type, investment type were also taken into account. Besides the above factor taken into consideration I have also taken into account the exporters relation with their buyer, input supplier, and inhouse process as other factors of competency. The analysis was done in the following way:   * Factor Analysis of factors retaining 1 to 10 scales   * Factor Analysis of factors retaining 1 to 5 scales   * Regression of FA of competition with the factors computed from the solution with 1 to 5 scales   * Factors of competition were multiplied with the exporters relation with buyer, input supplier and inhouse capacity, and then again regressed with the 1-5 factor solution.    * Competitive factors are taken as DV and factors 1 to 5 taken as IV Is there anybody who can tell me if I am following the right step for my analysis?


title: MCMC for Bayesian Inference (Gibbs Sampling) Varying Observed Data

text: At every step $k$, a Markov chain Monte Carlo algorithm for Bayesian inference with Gibbs sampling draws a parameter of the model to fit, $\beta_i^{(k)}$, from the conditional $Pr\left(\beta_i^{(k)}|\beta_1^{(k)},...,\beta_{i-1}^{(k)},\beta_{i+1}^{(k-1)},...,\beta_N^{(k-1)},X\right)$, where $X$ is the observed data. Do you know of some paper exploring this algorithm (or similia) where the data $X$ also vary along the chain? I do not mean at every step, just sometimes along the chain; imagine e.g. some financial transactions which happen almost continuously, and whose model always has to adapt. Thanks!


title: Convolutional networks: how to deal with duplicate filters?

text: I have to train a convolutional network. When I visualize its conv filter matrices from a layer, many filters are quite similar, almost equal. 1) Don't you know a common way to detect couples of similar kernels? Maybe, it is somehow connected with their correlation? 2) Don't you know any good strategies to deal with such duplicate kernels? My only offer is to re-initialize such duplicate with noise and continue the training, but I don't think it is the best idea.


title: How to compute asymptotic confidence intervals for differences in quantiles?

text: Can anyone give me advice on computing the asymptotic confidence intervals for a difference in quantiles of a distribution? For example, I have fit a log- normal distribution to doubly interval censored data with a maximum likelihood estimator and want to make statements about a percentage of the probability mass falling outside two fixed values (i.e. 23% (XX-YY) of values in this distribution fall outside of the range 1-2). I can get SEs for the median ($e^{\mu}$), dispersion ($e^s$), and quantiles ($e^{\mu+sq}$) using the delta method but for some reason can't think of how to get SEs for differences between them. To get the SEs for the quantiles and params above I use the following: $\bigtriangledown g(\mu, \log(\sigma)) = \left( \begin{array}{ccc} e^{\mu} & 0 & e^{\mu+s*z_{quants}}\\ 0 & e^{s+\log s} & z_{quants}e^{\mu +qs+\log s} \\ \end{array} \right)$ Which means that the new standard errors will be $\bigtriangledown g(\mu, \log(\sigma))^{T} H^{-1} \bigtriangledown g(\mu, \log(\sigma))$ where H is the approximated Hessian Matrix.


title: Power transformation on interactions in linear regression

text: In a multiple linear regression I have one predictor showing a cubic relationship with the outcome. If the same predictor $A$ is also part of an interaction with a continuous variable $B$, do I have to apply the polynomial transformation to both terms? $$ \hat{Y} = {c} + \beta_1{A} + \beta_2{A^2} + \beta_3{A^3} + \beta_4{A}{B} + \beta_5{A^2}{B} +\beta_6{A^3}{B}. $$ Logically it makes sense, but I've never seen a transformation applied to an interaction. Also I was wondering if the relationship can be simplified. I've found some information in this thread but my question is simply to understand if the transformation is correct.


title: How to show what is going on if we drop significiant variable in logit model?

text: How to show what is going on if we drop significant variable in logit model ? Bias and heteroskedasticity should emerge. But what is the framework for showing such behaviours in econometric models ?


title: Centering to use a sparse covariate?

text: In my dataset, there's a binary response, some factors, and some covariates. In particular, there are some covariates that are always present when `factor1=="A"` and these covariates are always missing (`NA`) when `factor1=="B"`. This missingness is structural, not random; it doesn't make sense for the covariate to have a value when `factor1=="B"`. To fit this as a logistic regression, I centered the covariates about the mean, and then set the covariate to 0 wherever it had been `NA`. I reason that in this way, the effect of the covariate on the response will be nullified whenever `factor1=="B"`, which is appropriate. Is this a reasonable approach to fitting these data? What is the name for this approach, and where (web, textbook) is there a good discussion?


title: Two samples, proportions

text: I have two groups of samples, each group of size 8. Each sample comes from an individual $i = 1 \dots 16$ and is a proportion $m_i$ of DNA methylated sites at a specific position in a genome. This can take a value between 0 and 1. We can safely assume that there is no variability between samples if they were taken from the same individual (i.e., the proportion in a single individual is calculated from a very large population of DNA strands). The individuals from both groups are matched. Now, more than one position has been tested (as above). Actually, almost half a million positions were tested, each of them corresponding to a proportion as described above. Effectively I have a matrix of the size $16 \times 5\cdot 10^5$, and each number in this matrix is a value between 0 and 1 corresponding to the proportion of DNA strands that are methylated in a specific sample at a specific genomic position. Given that I have a very large number of positions in a genome tested, I have an idea what the distribution of these measurements can be. Most of the sites are either methylated or not methylated -- the value $m_i$ is close to 1 or close to 0: ![Histogram of $m_i$ values](http://i.stack.imgur.com/0EfZz.png) H0 is that the groups do not differ. I have applied a Wilcoxon signed rank test and a t test (although clearly the assumptions of the t test do not hold, but it is generally considered robust), the results were similar. However, I wonder whether there might be better approach that takes the specific distribution of the $m_i$ values.


title: Error message from grm() in ltm: subscript out of bounds

text: I am analyzing a survey data set. One of the items in the survey is:               "I plan to work hard in my statistics course"       1=Strongly Disagree      2=Disagree      3=Neither Disagree or Agree      4=Agree      5=Strongly Agree      No one selected option 1 of Strongly Disagree. This is to be expected but the function `grm()` in package `ltm` generates an error: `"Subscript out of bounds"`. How do I get around this problem and still keeping the responses as they are or without removing that item?


title: What data structure is necessary for survival analysis?

text: I'm relatively new to survival analysis and try to get my data in the right shape. I have two tables both concerning the observed individuals. If I just would use one of the tables, I would have continuous information on each individual without any overlapping periods. As I however also need the information stored in the other table, it is necessary to merge the two tables. But then the episodes will be overlapping in some cases. I give you an example as illustration: Table 1:               ID: 1     start: 2000-01-01     end: 2002-12-31     state: A      Table 2:               ID: 1     start: 2002-01-01     end: 2002-04-15     state: Z      To do survival analysis (in Stata or R) does it matter if there are overlaps? If it does, do you have any suggestions on how to remove the overlaps?


title: Difference beween supplementary and active variables in PCA - Interpretation on obsevations?

text: I would like to introduce two supplementary variables into a PCA I'm conducting on a set of data measuring concentration in different material phases. However I'm unclear as to how to interpret the correlation of the "supplementary" variables to the principal component (as opposed to how I would interpret active variable loadings, and eventually the relationship between the PC and the scores). Clearly the factor loadings and sq. cosines of the variables are different depending on whether I treat them as active or supplementary. Specifically when variable "A" and "B" are included as active variables they have a strong correlation with PC2, but when treated as supplementary they do not link well to the first few PCs (and explain very little of the total variance). Based on my current understanding, the two different scenarios impact my interpretation of the scores. I realize this is likely a simple question, but would appreciate any help to clarify my understanding so I can make sure I proceed correctly. **UPDATE** To follow up then, one is investigating the relationships of the passive variables (eg: w1,w2,w3) with each other in the context of the PC structure. So if the loadings of all active variable (eg: v1, v2, v3) is high on PC1, but the passive variables have minimal calculated loadings on PC1 - how would one interpret the relationship between the passive variable and the active variable? There is no relationship between active and passive variables? **more info** here are poor quality charts - I will try and improve them. The first shows the red points as the active variables, and the second poor quality image shows the 3 variables treated as supplementary variables. I believe the scales are similar. The factor loadings and sq. cosines are different of course. ![active variables](http://i.stack.imgur.com/uvxZC.jpg) ![active and supplementary](http://i.stack.imgur.com/Mwi0E.jpg)


title: Do the predictions of a Random Forest model have a prediction interval?

text: If I run a `randomForest` model, I can then make predictions based on the model. Is there a way to get a prediction interval of each of the predictions such that I know how "sure" the model is of its answer. If this is possible is it simply based on the variability of the dependent variable for the whole model or will it have wider and narrower intervals depending on the particular decision tree that was followed for a particular prediction?


title: Clopper-Pearson for non mathematicians

text: I was wondering if anyone can explain to me the intuition beyond the Clopper- Pearson CI for proportions. As far as I know, every CI includes a variance in it. However, for proportions, even if my proportion is 0 or 1 (0% or 100%), the Clopper-Pearson CI can be calculated. I tried looking at the formulas, and I understand it has something with percentiles of the Binomial distribution and I understand that finding the CI involves iterations, but I wondered if anyone can explain the logic and rational in "simple words", or with minimum math ?


title: Random Forests and data transformations

text: I've started reading about Random Forests and one of the attributes that appeals to me is that they are good at dealing with independent variables that interact with one another. Does Random Forest also "automatically" transform variables into, for example, the square of the variable? Is there a good reason to square or perform other transformations of a variable (of course dependent on the situation)? For example if I have `y = x1 + x2 + x3` is there anything to be gained by running the model `y = x1 + x1^2 + x2 + x2^2 + x3 + x3^2`?


title: chi-squared test for non-iid binomial data

text: **Simple case:** Suppose we have $N$ observations of $X_i$ that we believe are theoretically iid binomial random variable $\text{Binomial}(k, p)$. Well, one way we could test to see if this is true is using the standard Pearson's chi- squared statistic. **The question:** Suppose we have a more difficult problem, where we have $N$ different $X_i$'s, but now each is theoretically from $\text{Binomial}(k, p_i)$, with constant $k$ between all draws, _and_ we know each $p_i$ corresponding to each observation. Is there an obvious analytic extension of the previous test to this case? At the moment I am simulating a null distribution of chi-squared statistics from the known $p_i$'s and testing from that.


title: Why do coefficients and significance levels change so much in my OLS? Do I need MIXED model instead?

text: Please consider the following OLS model: ![enter image description here](http://i.stack.imgur.com/9l9iB.png) As a rule of thumb, a VIF score of max. 5 is considered acceptable collinearity in most fields. So for the better part I believed this model to be quite nice and significant (R-squared .516). However, by adding only 2 variables (at the bottom of the model), many values change significantly. For instance, `GCI-505` goes from being significant at the 10% level (p=.008) to being insignificant (p=.677). `INCID` reacts similarly, whereas `HOF_PDI` and `HOF_LTO` have a reverse change: ![enter image description here](http://i.stack.imgur.com/3jy4O.png) My question is: what does this mean for my model? It doesn't seem very robust, although I though that it should have been with a significant F-score, and not too high VIF scores. What am I missing / lacking in knowledge? `[edit 1]` Some explaining of the variables: the model tries to explain the extent to which a company is behaving socially responsible. The dataset holds 4,000 companies from 35 countries. First, each company has a CSR score (the DV), a marketvalue and bookvalue (`MKT_AVG_LN` and `PTBV_AVG_LN)`, a number of treaties on responsible investment they have signed (`SIGN`), and a number of incidents that have been reported in the media (`INCID`). Second, each company is headquartered in one of 35 countries. These countries all have certain properties, for instance from which legal family their legal system stems and the average level of education. These country-specific values have been included as a variable for each company. For instance, a company with `LAW_FRA`-dummy set to 1 is headquartered in a country with a french legal system. When a company has `LAW_SR` set to `5.45`, this means that the shareholder rights index for the country the company in headquartered in, is `5.45`. `LAW_FRA`, `LAW_SCA` and `LAW_ENG` (not included in the model) are dummies for legal family. The rest of the variables are continuous. There are possibly 2 problems with my dataset   * I have "elevated" country-variables to a company level;   * the number of companies from each country is very unequally divided (for instance, 40% comes from the USA whereas only 1% of companies is headquartered in Singapore. This means that if USA has a `LAW_SR` score of `5.45`, 40% of companies in the dataset have a `LAW_SR` score of `5.45` whereas only a few companies will have a `LAW_SR` value of `3.00` (the `LAW_SR` score for Singapore). `[edit 2]` In SPSS Mixed Models has the following dialogue: ![enter image description here](http://i.stack.imgur.com/dnjjI.png) I am confused about which variables are factor, covariate and residual weight. As I believe continuous company-variables, such as a companies' market value and number of signed treaties (`MKT_AVG_LN` and `SIGN`) are factors, and company variables that are lifted from a country-level (for instance, `LAW_SR`, which is a numerical value of the relative strength of shareholders under company law of a country) are covariates. Is this correct? Also, I have a few dummy variables. For instance, each company operates under either French, German, English or Scandinavian law. Here, `LAW_FRA`, `LAW_GER`, `LAW_ENG`, `LAW_SCA` are dummy variables indicating for each company the origin of company law in the country where the company is headquartered. These dummy values are derived from the variable `LAW_FAMILY`, which for each company is set to either `FRA`, `GER`, `ENG` or `SCA`. In mixed models, should I discard the dummy variables and use `LAW_FAMILY` as a covariate also?


title: Comparing the magnitude of the standard deviation with the average value

text: For calculating the concentration of an unknown sample, my average value after three trials was about .03425 +- twice the standard deviation (.00159) In comparing twice the standard deviation to the average value, is my precision about average?


title: Use the Central Limit Theorem to calculate the approx probabilities of Gamma RVs?

text: If $X_i, i=1,...,$ are independent, identically distributed $\operatorname{N}(0,1)$ random variables, and $Y_i = X_i^2$ are independent $\operatorname{Gamma}(\frac{1}{2},\frac{1}{2})$ RVs, use the Central Limit Theorem to calculate the approximate values of the probabilities:      $\sum_{1}^{50} {Y_i \gt 40}$ and $\sum_{1}^{50} {Y_i \lt 60}$      I'm not really sure how to go about this.


title: How do I vertically stack two graphs with the same x scale, but a different y scale in R?

text: Greetings, Currently I'm doing the following in R:               require(zoo)     data <- read.csv(file="summary.csv",sep=",",head=TRUE)     cum  = zoo(data$dcomp, as.Date(data$date))     data = zoo(data$compressed, as.Date(data$date))     data <- aggregate(data, identity, tail, 1)     cum  <- aggregate(cum, identity, sum, 1)     days = seq(start(data), end(data), "day")     data2 = na.locf(merge(data, zoo(,days)))          plot(data2,xlab='',ylab='compressed bytes',col=rgb(0.18,0.34,0.55))     lines(cum,type="h",col=rgb(0,0.5,0))      Snip of summary.csv:               date,revision,file,lines,nclass,nattr,nrel,bytes,compressed,diff,dcomp     2007-07-25,16,model.xml,96,11,22,5,4035,991,0,0     2007-07-27,17,model.xml,115,16,26,6,4740,1056,53,777     2007-08-09,18,model.xml,106,16,26,7,4966,1136,47,761     2007-08-10,19,model.xml,106,16,26,7,4968,1150,4,202     2007-09-06,81,model.xml,111,16,26,7,5110,1167,13,258     ...      The last two lines plot the information I need, and the result resembles the following: ![alt text](http://i.stack.imgur.com/ech9l.png) Blue line is the entropy in bytes of the artifact I'm interested. Green lines represent the entropy of the changes. Now, in this graph, it works well because there isn't a huge difference in scales. But I have other graphs where the green lines become so small one cannot see. The solution I was looking for, involved two things:   1. To move the green vertical lines to a second graph, just below the first one, with its own y axis, but shared x axis.   2. To provide it a logarithmic scale, since I'm more interested in the "magnitude", than in the specific values. Thanks in advance! P.S. If someone can also tell me how could I put "minor ticks" in the x scale referring to the months, I appreciate :-) If these are too much questions for a single post, I can divide them further.


title: Speed up web a/b tests with sample size checkpoints

text: Before starting an a/b test with a control and one experimental route, I can calculate a required sample size based on conversion rate estimates for both routes. I can get a good estimate of the conversion rate of the control by looking at historical data. But the conversion rate of the experimental route is unknown. What I would like to do is calculate a number of different sample sizes based on a variety of sensitivities. For example, I can calculate sample sizes for a 10%, 15% and 20% sensitivity (increase in conversion from the control) that might look like this:               Sensitivity   Required Sample Size     10%           1,961     15%           871     20%           490      Some of the reading I've done says that you should calculate a single sample size at the start of the test and always run the test for that long. Question:   * Is there any problem with checking for statistical significance at multiple pre-calculated sample sizes and potentially ending a test early if I've found the results to be statistically significant?  Example: I originally estimate that the experimental route will outperform the control by 15%. But once I've reached 490 samples I find that the experimental route is actually outperforming the control by 20%, can I end the test and declare that the experimental route boosts conversion by 20%?


title: What are the independence assumptions of Cox models with time varying covariates?

text: In longitudinal studies, you might be observing time-to-event endpoints with some covariate that changes as a function of time. When covariates are fixed at baseline, the only independence assumption involved is within the risk set, or the collective at risk observed during event times. Otherwise, the fact that an individual survived at time 4 is in no way expected to influence whether that individual survived at time 8 (aside from the modeled covariates and baseline hazard). When individuals have a discrete change in covariates, they are censored at that time and reintroduced into the analysis at another point. In this case, I expect there should be no reason to handle the fact that individuals were in risk sets at another point in time with different values. The only exception I can conceive of is whether the influence of the covariate has any cumulative effect (but that seems to be a significantly different question, and more of an assumption about the linear model than Cox regression). Can you confirm that my reasoning is correct?


title: Given this time series, what statistical methods would be used for description and forecasting?

text: These static cumulative default rate tables and charts come from this public report published by a credit rating agency. Basically, you take all the loans originated in a period of time (a "vintage") and track the cumulative amount of money due but not repayed by the borrowers as a percentage of the amount of money lent (the "cumulative default rate") for each period of time. Each curve in the chart (each vintage) shows the behaviour of _all_ the credits given by the bank during a given period of time (say, a month or a quarter). This mean:   * Observations included in all the curves are originated by the same bank using the same rules, therefore they share common characteristics and similar behaviour.   * The dataset analyzed includes all the loans originated by the bank. The observations include the whole population, they do not come from a sampling process. Some visualization: **The dataset**               Vintages Amount 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21      1Q01 612,000 0.10 0.57 1.13 1.89 2.74 3.34 3.61 4.06 4.36 4.66 4.73 4.83 5.08 5.18 5.27 5.38 5.39 5.41 5.43 5.48 5.58      2Q01 700,000 0.09 0.64 0.99 1.77 2.58 3.03 3.46 3.64 4.01 4.28 4.32 4.47 4.73 4.82 4.97 5.05 5.07 5.10 5.23 5.32      3Q01 732,490 0.13 0.70 1.53 2.08 2.60 3.08 3.28 3.75 3.82 4.32 4.37 4.57 4.70 4.92 4.98 4.99 5.07 5.12 5.14      4Q01 793,984 0.23 0.85 1.73 2.19 2.72 3.31 3.79 3.96 4.40 4.66 4.78 4.92 5.13 5.21 5.37 5.38 5.49 5.57      1Q02 834,563 0.22 0.73 1.28 1.89 2.77 3.11 3.45 3.82 3.91 4.37 4.45 4.63 4.85 5.08 5.11 5.15 5.16      2Q02 904,485 0.14 0.65 1.45 2.17 2.72 3.20 3.81 4.08 4.38 4.43 4.54 4.91 5.08 5.11 5.15 5.18      3Q02 953,534 0.10 0.71 1.42 1.98 2.55 3.14 3.71 4.15 4.42 4.67 4.73 4.89 4.95 5.03 5.24      4Q02 980,124 0.13 0.64 1.13 1.61 1.98 2.23 2.71 3.03 3.30 3.68 3.77 3.82 3.87 4.08      1Q03 1,202,324 0.18 0.69 1.14 1.95 2.42 2.56 2.98 3.15 3.67 3.68 3.76 4.00 4.10      2Q03 1,014,485 0.16 0.59 1.39 1.90 2.23 2.58 2.92 3.46 3.68 3.74 3.87 3.99      3Q03 1,171,731 0.15 0.60 1.18 1.66 1.97 2.45 3.11 3.15 3.41 3.50 3.81      4Q03 1,312,338 0.13 0.61 1.13 1.44 2.11 2.33 2.85 2.87 3.01 3.05      1Q04 1,509,189 0.16 0.59 1.14 1.45 1.75 1.99 2.36 2.47 2.63      2Q04 1,569,557 0.13 0.58 1.09 1.46 1.85 2.19 2.28 2.44      3Q04 1,545,922 0.14 0.58 0.99 1.39 1.72 1.97 2.18      4Q04 1,514,485 0.18 0.63 0.94 1.46 1.81 1.84      1Q05 1,574,685 0.16 0.55 1.04 1.45 1.72      2Q05 1,578,485 0.10 0.38 0.74 0.90      3Q05 1,501,234 0.06 0.44 0.64      4Q05 1,472,187 0.07 0.41      1Q06 1,459,861 0.08       ![enter image description here](http://i.stack.imgur.com/SrTIM.png) **The chart** ![enter image description here](http://i.stack.imgur.com/r3Erw.png) **For better visualization, I add a chart from another dataset with more observations.** ![enter image description here](http://i.stack.imgur.com/EGwpB.png) Since I'm not an expert in statistics, I'd like to know what kind of statistical methods I should study to be able to a) describe the behaviour of each vintage (for example, indicators to compare new credits versus old credits), and b) forecast the final cumulative default rate (the default rate at the end of the life of the credits) for the youngest vintages. I know the I've only included very little information about the case, but I'm willing to do my research on any statistical method you would suggest as possibly useful.


title: Survey questionnaire rebasing after removal of a selection

text: Suppose that I have a market research survey with the question "What brand of television are you planning on buying" and then also have a choice that is "I don't plan on buying a television." Respondents may choose more than one brand, but they may not make any other choices if they select that they do not plan on buying a television. **If I wanted to do an analysis of television brands among only those people who said they were going to buy a television, can I remove all of the respondents that said they would not buy a television, then recalculate the proportions for each brand with the new base number of observations?** If this skews the demographic distribution of my sample, is that okay? (e.g. the full sample has 30% Caucasians and 10% Asians, but the sample of people who say they will buy televisions is 30% Asians and 10% Caucasians. It seems like this should be okay, since we would not expect the two populations (TV buyers and people in general) to be the same.


title: What are the pitfalls of using subjective ratings for correlation and causal explanation?

text: Take a study that collects subjective guesstimates as a proxy for some variable, like the quality of service. In a simplest model, we then regress the monthly sales on these guesstimates. How to interpret the results given that they originate from subjective estimations? One problem is that respondents have no common anchor in many cases (as explained on CogSci.SE). This means that even if the relations between the quality of service and sales exist, they may not appear because of difference in understanding what "good" service is. Let's take the extreme case when the sales and quality are the same across the sample, but each person rates the quality as he or she perceives it. Then we have a non-existing variation in the quality. What are the techniques to normalize subjective scales? Kahneman and Tversky suggested the reference class model, in which you adjust the cognitive bias using a common class for calibration across the whole sample of experts. But it's not always feasible.


title: Calculate coefficients in a logistic regression with R

text: In a multiple linear regression it is possible to find out the coeffient with the following formula. $b = (X'X)^{-1}(X')Y$               beta = solve(t(X) %*% X) %*% (t(X) %*% Y) ; beta      For instance:               > y <- c(9.3, 4.8, 8.9, 6.5, 4.2, 6.2, 7.4, 6, 7.6, 6.1)     > x0 <- c(1,1,1,1,1,1,1,1,1,1)      > x1 <-  c(100,50,100,100,50,80,75,65,90,90)     > x2 <- c(4,3,4,2,2,2,3,4,3,2)     > Y <- as.matrix(y)     > X <- as.matrix(cbind(x0,x1,x2))          > beta = solve(t(X) %*% X) %*% (t(X) %*% Y);beta              [,1]     x0 -0.8687015     x1  0.0611346     x2  0.9234254     > model <- lm(y~+x1+x2) ; model$coefficients     (Intercept)          x1          x2       -0.8687015   0.0611346   0.9234254       I would like how to calculate in the same "manual" way the beta for a logistic regression. Where of course the y would be 1 or 0. Assuming I'm using the binomial family with a logit link.


title: Is it appropriate to treat n-point Likert scale data as n trials from a binomial process?

text: I've never liked how people typically analyze data from Likert scales as if error were continuous & Gaussian when there are reasonable expectations that these assumptions are violated at least at the extremes of the scales. What do you think of the following alternative: If the response takes value $k$ on an $n$-point scale, expand that data to $n$ trials, $k$ of which have the value 1 and $n-k$ of which have the value 0. Thus, we're treating response on a Likert scale as if it is the overt aggregate of a covert series of binomial trials (in fact, from a cognitive science perspective, this is actually an appealing model for the mechanisms involved in such decision making scenarios). With the expanded data, you can now use a mixed effects model specifying respondent as a random effect (also question as a random effect if you have multiple questions) and using the binomial link function to specify the error distribution. Can anyone see any assumption violations or other detrimental aspects of this approach?


title: What's some examples of the fact disjoint events are dependent?

text: Disjoint events A and B with positive probability must not be independent as $P(A)$, $P(B) \gt 0$, $P(A\cap B) = P(\emptyset) = 0 \Rightarrow P(A\cap B) \lt P(A)P(B)$. What are some real world examples that can help better understanding this?


title: Interpretation of survival risk predictions generated from R's predictSurvProb

text: I am starting to learn R and its power for giving survival predictions. ### Abstract I use the `predictSurvProb` function from the `pec` package to get predicted survival risks. My data has 5000 rows wherein each thousand falls into a certain `time` category. These `time` categories are separated by two-year intervals; thus for time periods one to five, these have 1000 records that fall within these ranges. My data has 47 variables: the time and event; 20 numeric, 10 binary and 15 nominal predictive variables. My steps to generate the risks are as follows:   1. Fit a Cox model, `xn`s are the rest of the predictive variables.              coxmodel <- cph(Surv(time,event)~x1+x2+x3+x4+x5,          data=myData, surv=TRUE, x=TRUE, y=TRUE)         2. Declare time periods.              myTime <- c(24,48,72,96,120) # months     periods <- quantile(myTime, conf.int=TRUE)        3. Get the prediction risks!              survRiskPreds <- predictSurvProb(coxmodel, newdata=myData, times=periods)     round(preds, digits=6)      ### Problem (Note that these below are example results but nonetheless represent the phenomena that I'm facing.) Since I have five different time periods, I expect each of the 1000 records to have optimistic results in their own time period for example:                                                 24      48      72      96     120      Record survives for 30 months : 0.96251 0.81881 0.65141 0.46836 0.11000     Record survives for 60 months : ?       0.96251 0.81881 0.65141 0.46836     Record survives for 90 months : ?       ?       0.96251 0.81881 0.65141     Record survives for 100 months: ?       ?       ?       0.96251 0.81881      Record survives for 150 months: ?       ?       ?       ?       0.96251      But in reality the results have the pattern for all records no matter what their `time` value is:                                                 24      48      72      96     120      Record survives for 30 months : 0.96251 0.81881 0.65141 0.46836 0.11000     Record survives for 60 months : 0.98111 0.83809 0.71174 0.60173 0.21272     Record survives for 90 months : 0.95411 0.82209 0.69974 0.51936 0.20272     Record survives for 100 months: 0.99911 0.88909 0.60074 0.57736 0.21272     Record survives for 150 months: 0.98801 0.78109 0.68211 0.51436 0.33272      ### Sample data and code I made a dummy dataset that more or less represents the one I currently have. Below is the script I use to manipulate the data. I included all predictive variables in the Cox model.               library(survival)     library(pec)      library(rms)          theData <- read.csv("dummy.csv",      colClasses=c(        "female"="factor",        "bin1"="factor",        "bin2"="factor",        "marital"="factor",        "yrDiagnosis"="factor",        "response1"="factor",        "response2"="factor",        "response3"="factor",        "country"="factor",        "race"="factor"))     attach(theData)          # removed one predictive variable (yrDiagnosis) because it caused an error     coxmodel <- cph(Surv(time,event)~age+num1+num2+num3+num4+       num5+female+bin1+bin2+marital+response1+       response2+response3+response4+response5+response6+country+race         , data=myData, surv=TRUE, x=TRUE, y=TRUE)      coxmodel           myTime <- c(24,48,72,96,120)     periods <- quantile(myTime, conf.int=TRUE)     periods          preds <- predictSurvProb(coxmodel, newdata=theData, times=periods)     round(preds[1:5,], digits=6)       # survival within 24 months     round(preds[1000:1005,], digits=6) # survival within 48 months     round(preds[2000:2005,], digits=6) # survival within 72 months     round(preds[3000:3005,], digits=6) # survival within 96 months     round(preds[4000:4005,], digits=6) # survival within 120 months      * * * Am I doing something wrong ie. wrong format in the `predictSurvProb` steps? I am not familiar with other R packages but I would appreciate if anyone guides me get to my desired result.


title: Fitting probability distribution to data

text: ![enter image description here](http://i.stack.imgur.com/0jywZ.png) I am trying to fit a model for the values plotted above. The explanatory variable represents amounts of compensation claim in an earthquake, and the response variable represents amounts of compensation awarded. Can someone tell me what probability distributions would be my possible options for these data? Thanks in advance.


title: Calculating Jaccard or other association coefficient for binary data using matrix multiplication

text: I want to know is their any possible way to calculate Jaccard coefficient using matrix multiplication. I used this code                   jaccard_sim <- function(x) {         # initialize similarity matrix         m <- matrix(NA, nrow=ncol(x),ncol=ncol(x),dimnames=list(colnames(x),colnames(x)))         jaccard <- as.data.frame(m)              for(i in 1:ncol(x)) {          for(j in i:ncol(x)) {             jaccard[i,j]= length(which(x[,i] & x[,j])) / length(which(x[,i] | x[,j]))             jaccard[j,i]=jaccard[i,j]                    }          }      This is quite ok to implement in R .I have done the dice similarity one, but got stuck with tanimoto/jaccard. Anybody can help?


title: Errors fitting a Penalized Linear Model in R

text: I am trying to reproduce the results from this answer in CrossValidated with no luck. I keep getting an error when calling:               print("Fitting Penalized Linear Model")     GLMmodel <- train(X,Y,method='glmnet',trControl=myCtrl,tuneLength=10)      Here is the error:               [1] "Fitting Penalized Linear Model"     Error in { :        task 1 failed - "arguments imply differing number of rows: 1264, 736"     In addition: There were 50 or more warnings (use warnings() to see the first      50)                                                                             The full code is below:               #Setup     rm(list = ls(all = TRUE))     set.seed(1)          #Generate an input matrix with 10,000 values     #2,000 rows, 5 columns     X <- data.frame(replicate(5, rnorm(2000)))          #Construct Y using X     Y <- runif(1)*X[,1]*X[,2]^2+runif(1)*X[,3]/X[,4]          #Convert Y to binary     Y <- as.factor(ifelse(sign(Y)>0,'X1','X0'))          #Create bootstrap samples for fitting models     library(caret)     print("Creating bootstrap samples")     tmp <- createResample(Y,times = 25)     myCtrl <- trainControl(method = "boot", index = tmp, timingSamps = 10)          #Fit models     print("Fitting Random Forests")     RFmodel <- train(X,Y,method='rf',trControl=myCtrl,tuneLength=1)                      print("Fitting SVM")     SVMmodel <- train(X,Y,method='svmRadial',trControl=myCtrl,tuneLength=3)          print("Fitting Penalized Linear Model")     GLMmodel <- train(X,Y,method='glmnet',trControl=myCtrl,tuneLength=10)          print("Fitting Nearest neighbors")     KNNmodel <- train(X,Y,method='knn',trControl=myCtrl,tuneLength=10)          print("Fitting Neural Network")     NNmodel <- train(X,Y,method='nnet',trControl=myCtrl,tuneLength=3, trace = FALSE)          #Assess re-sampled (out of sample) accuracy     print("Assessing re-sampled (OOB) accuracy")     resamps <- resamples(             list(   RF = RFmodel,                     SVM = SVMmodel,                     GLMnet = GLMmodel,                     KNN = KNNmodel,                     NN = NNmodel                     ))     print("Summarizing resamps")     summary(resamps)          print("Plotting accuracy")     dotplot(resamps, metric = "Accuracy")          #Assess run-time (in seconds)     print("Assessing timings")     resamps$timings      I am relatively new to R so I can't tell what's wrong with that line. Did the package GLM change since the time of that post? Is this a possible bug with `GML`? For reference. I am running the latest version of R `(2.15.2)` with the latest stable version of caret `(5.15-052)`


title: Self-Study question on Normal Distribution

text: I was attempting a self-study question on the normal distribution, but I wasn't able to get to the answer. I'm not sure if my method used here is incorrect, and would appreciate some guidance and advice. **Question** > The pulse rate per minute of the adult male population between 18 and 25 > years of age in the United States is known to have a normal distribution > with a mean of 72 beats per minute and a standard deviation of 9.62. The > requirements for military service state that anyone with a pulse rate over > 100 is medically unsuitable for service. In a group of 20 males between 18 > and 25 years of age, what is the probability that at least one of them would > be declared unfit because their pulse rates are too high? Assume that their > pulse rates are independent. **My approach was as follows:** > $X \sim \mathcal{N}(72, 9.62^2)$ > > With $n = 20$, find $\Pr(X>100)$. Since the value is normally distributed, I > choose to normalize the information and use $Z$ to find my value. > > $\Pr[Z > (100-72)/9.62] = \Pr(Z > 2.9106)$. Using the $Z$-table, I deduced > that the value is 0.00181. **_This answer, though, is far off from in the guide's answer key, which is 0.0356. The answer key does not work through the steps, so there is no information about where I went wrong._** However, this answer seemed to be wrong based on the study guide. Appreciate some guidance and directions please.


title: Question on Probability

text: I was attempting to solve a question on finding the probability that a randomly selected person had an income of $50,000 or more would be audited, based on information on the table below: ![enter image description here](http://i.stack.imgur.com/0aTNw.jpg) My attempt was to take the proportion of those with income more than $50,000 and multiply it with the percentage audited. > **Proportion of those earning more than $50,000** : 11.4 / 101.5 = 0.11231 > > 0.11231 * 2.24 = 0.2516 However, my answer was slighty off from the provided answer of 0.2394. Am I even on the right track with the method that I am using? Appreciate any advice please.


title: Expected value of item in a sorted list of integers

text: If you were to take random positive integers, put them into a list, and sort them is there anyway to find the expected value of the kth item in the list? The list is sorted in ascending order. By random positive integer I mean some number in the range 0-N where the probability of each number occurring is equal. (I think what I'm trying to say is uniform distribution).


title: Difference between regression methods

text: I have a set of data where the response is a proportion. For each event in the experiment, a system will correctly tag `X` of `Y` items. In the end, I will have something like 100 different results. I have been trying to determine what is the most correct method of fitting the response variable. So far I have uncovered 4 number of options.   1. R's built-in `glm` using one line for each of the `Y` items in each experimental event   2. R's built-in `glm` using a proportion (`X/Y`) with the weight set to `Y` for each experimental event   3. The `betareg` fit method from the R `betareg` package   4. Quasi-likliehood method proposed by Papke and Wooldridge (pdf) * * * I am unsure how to rank these different methods to know which is the most appropriate. Here are some specific questions:   1. How does the likelihood change for these different options (especially between options 1 and 2)?   2. What are the pros and cons of the different methods?   3. Do these methods have limitations on the types of factors that can be built into the model?   4. For the `betareg` function, are the weights specified the same way as the `glm` function?   5. Most people would tell me to use option 2 listed above, but I have seen a lot of pull for option 3. Why should I choose option 3 over option 2?   6. What do you do if you do not know the denominator of the proportion (i.e., you cannot properly weight the data for option 2)? As you can see these questions are quite general and could be the elements of a full academic text. I was just not sure where to go to get the answers to these questions.