title: Design matrix versus contrast matrix

text: Suppose X is the design matrix of my experiment of which I want to model a linear regression model. Such a design matrix can be created by (in this example it is a full factorial design)               library(BHH2)     Des2 <- ffDesMatrix(5)      Now, a contrast matrix can be obtained by the following code, corresponding to ( t(X) %x% X)^-1 %x% t(X) (in which %x% represents matrix multiplication):               solve(t(Des2) %*% Des2) %*% t(Des2)      This results in a matrix giving the Y observations that will be contrasted against each other in order to obtain a parameter estimate (e.g. first row will correspond to beta1). Now, one can also make a design matrix using the model.matrix function               Des <- as.data.frame(Des2)     names(Des) <- c("X1","X2","X3","X4","X5")     model.matrix(~X1*X2*X3*X4*X5,data=Des)      This is a model matrix that can be used to estimate all possible interactions between five factorial variables with two levels. When comparing the columns of the model matrix obtained by the last piece of code with the rows of our contrast matrix, one can note that they are rather similar. The values might be different, but the signs representing the contrasts are identical for a row/column representing the same parameter to be estimated. My question now refers to these columns and rows: is there any difference to be noted between these two? Or are they identical as in that they both represent contrasts of the Y observations that will be used for estimation of the parameter, and have I thus summed up two equivalent way in obtaining these contrasts?


title: Relation between Gaussian mixture models and maximum likelihood?

text: I need some help understanding the relation between the maximum likelihood and Gaussian mixture models. I have seen that there is a relationship between the expectation maximization algorithms and these two terms, but I could not quite understand the relation well enough. Could some please give me an example of how these two terms are related or different?


title: Long tail of pareto vs. pearson IV distribution

text: I am told that wealth follows a Pareto distribution, and that IQ follows a Pearson IV distribution (http://www.abelard.org/burt/burt-ie.asp). Both Pareto and Pearson IV distributions have long tails. Which one is longer? In plain language: are billionaires further from the average population in terms of wealth than geniuses are in terms of IQ? And considering that wealth is a cardinal variable and IQ and ordinal one, does this question even make sense?


title: OLS Coefficient estimator; Transformation from Matrix to sum of matrices form

text: I do not understand why the following equality holds (taken from Cameron & Trivedi 2005: _Microeconomtrics_ ): $\hat{\underline{\beta}}_{OLS}=(\textbf{X}'\textbf{X})^{-1}\textbf{X}'\textbf{y}=(\sum_{i=1}^{n}\textbf{x}_i\textbf{x}_{i}^{'})^{-1}\sum_{i=1}^{n}\textbf{x}_iy_i$ The notation is $ \textbf{x}_{i}^{'} =(x_{i1},...,x_{ik})$ where $i$ designates the different observations of $k$ variables and $\textbf{X}$ stacks these upon another. I have tried to develop the RHS of the equality, but I keep ending up with a scalar instead of a vector of coefficients. Any help would be greatly appreciated!


title: How do you suppress imputed data and only display results for the pooled data?

text: I've done 15 imputations and am now running various tests on the new data, but all I'm really interested in is the pooled results at the end. The long list of imputations (especially in a test like an independent T) is unnecessary and I'd like to just hide all of it/suppress it when I run any tests on imputed data. is there anyway to do that in SPSS Statistics 22?


title: Estimating the effect of one time series on another in the context of personal health

text: **Short version:** I have two time series: steps taken and change in body fat mass. Both are daily data. I am trying the estimate the effect of the number of steps taken on the change in body fat mass. I am especially interested in getting a point estimate and interval for for where the effect goes from positive to negative. **Loong version (including background, details, and methodological ramblings):** A few months ago I was reading a bunch of health related blogs and came across some references to the exercise physiology literature indicating that it is not the specific amount of exercise that has an effect on obesity, but rather the general level of sedentary time of the person. What matters is the integral of activity over time, and exercising a couple of times a week does not have a significant effect on it. I am the kind of person who does essentially no exercise at all, but I do tend to get restless and move around when sitting in place long stretches of time. So I decided to experiment on myself to quantify how my amount of sedentary time impacts my body fat. The problem with relying strictly on body mass is that there are other things that contribute to it, like water weight or muscle mass, that could reasonably be expected to change depending on the amount of sedentary time. My solution to this was to buy a scale that measures both weight and body fat percentage. By multiplying weight with body mass percentage (divided by 100) I get body fat mass. The proper way to measure sedentary time is with a three axis accelerometer, but a good one is over $500 and so was outside my budget. I had to make do with a high end step counter. The problem with a step counter is that it doesn't differentiate between running fast up a hill and gently strolling down a hill. But I am not the kind of person to be running up hills, nor do I live in a hilly area, so I judged it shouldn't be a problem. So I set about measuring these things. Additionally I set myself the goal to walk places, to that I get at least one hour of walking time per day. But being a lazy person weak of will I would only achieve that goal intermittently. That seemed to me a good thing in that I would get a lot of variability in the number of steps taken to work with. I have been measuring my weight and body fat percentage every morning after hitting the toilet (to minimize the effect of the weight of waste in the body) and measured the number of steps taken for a few months now. Being bot curious and weak willed of nature I wanted to get some preliminary analysis done now. The question then becomes how to analyze this data I have so far collected, which, after filtering out missing data points, adds up to 87 data points. I thought about simply running a simple regression and then solving for zero in the resulting model. That would give a a point estimate of where the number of steps start having a negative effect on body fat mass. But I am unsure if the time-series nature of the data would invalidate the result. Nor do I have any idea how I would get a confidence interval for that estimate. I realize I've rambled on quite a bit here but I am not exactly sure which details are relevant and should be taken into account in the analysis. **Edit:** steps:               8999, 10823, 10025, 4282, 7072, 5895, 8240, 7875, 14176,      9512, 9454, 5854, 1648, 1834, 10291, 2368, 2884, 7767,      8026, 1742, 1745, 2629, 8452, 6067, 6215, 6502, 10367,      7464, 4120, 9644, 5684, 8990, 5446, 8777, 8799, 8100, 8904,      4846, 4283, 7276, 1784, 6343, 7635, 12544, 3644, 3340,      4244, 12060, 5485, 6928, 3158, 9358, 5015, 10077, 7988,      8329, 5954, 2237, 4753, 5992, 6982, 7527, 8813, 4438, 8426,      7926, 6465, 7660, 8254, 7354, 1032, 6417, 4939, 7562, 8789,      3895, 3273, 3364, 4358, 8873, 8512, 7248, 4215, 1058, 3904,      8309, 7159      body fat mass:               17.38, 17.3383, 16.0398, 16.758, 15.9996, 16.1398, 16.8328, 16.3385,      15.1272, 16.2588, 16.6014, 16.5776, 16.4565, 16.2155, 16.4979,      15.9984, 16.4358, 16.5528, 15.9594, 15.6418, 15.9594, 16.3812,      17.1785, 16.863, 16.744, 16.9812, 16.842, 16.6221, 17.114, 15.9001,      16.359, 16.9812, 17.2584, 17.4618, 16.4016, 16.744, 17.4528, 16.5444,      17.1842, 17.0826, 17.4618, 16.3358, 15.8004, 16.5737, 17.3383,      16.4736, 16.1588, 17.114, 15.6618, 16.8167, 16.6155, 16.9974,      17.1039, 16.9916, 17.91, 16.9644, 15.9399, 17.3528, 17.064, 17.2791,      15.678, 16.1385, 16.5946, 16.116, 16.611, 16.6779, 17.292, 16.116,      14.0504, 16.4151, 16.5737, 17.446, 14.6306, 16.4151, 16.6782, 16.464,      16.669, 16.5424, 16.2288, 16.277, 16.695, 16.277, 16.7323, 17.1864,      17.1072, 16.653, 15.9996      At least where body fat mass is concerned I figured the mass one day would essentially be the mass of the previous day plus some small change, so it might make more sense to take the differences of the mass and model using that.


title: Analysis of Split (A/B) tests using Poisson and/or Binomial Distribution

text: I've seen this question asked before, but I have yet to come across a definitive source answering the specific questions:   * What's the most appropriate statistical test to apply to a small A/B test?   * What's the R code and interpretation to analyze a small A/B test? I'm running a small test to figure out which ads perform better. I have the following results: **Position 1:** `variation,impressions,clicks row-1,753,26 row-3,767 7 ` **Position 2:** `variation,impressions,clicks row-1,753,16 row-3,767 13 ` **Position 3:** `variation,impressions,clicks row-1,753,2 row-3,767 7 ` I think it's safe to say these numbers are small and likely to be not normally distributed. Also, it's click data so there's a binary outcome of clicked or not and the trials are independent. **Appropriate test** In analyzing each position for significance, I think comparison with a binomial or Poisson distribution makes the most sense. According to the OpenIntro Stats (and other sources) book, a variable follows a Poisson distribution "... if the event being considered is rare, the population is large, and the events occur independently of each other." The same source classifies a binomial variable approximately the same way adding that the probability of success is the same and the number of trials is fixed. I appreciate this is not an either/or decision and analysis can be done using both distributions. Given A/B (split) testing is a science that has been practiced for several years, I imagine that there is a canonical test. However, looking around the internet, I mostly come across analysis that uses the standard normal distribution. That just seems wrong :) Is there a canonical test to use for A/B tests with small #'s of clicks? **Interpretation and R code** I've used the following R code to test significance for each position: Position 1:               binom.test(7, 767, p=(26/753))          Exact binomial test          data:  7 and 767     number of successes = 7, number of trials = 767, p-value = 1.077e-05     alternative hypothesis: true probability of success is not equal to 0.03452855     95 percent confidence interval:      0.003676962 0.018713125     sample estimates:     probability of success                 0.009126467      I interpret this result to mean: The probability of success in the test group is indeed different than the control group with a 95% confidence interval that the success probability is between .368% and 1.87%               ppois(((26-1)/753), lambda=(7/767), lower.tail = F)     [1] 0.009084947      I interpret this result to mean: Given a Poisson distribution with a click rate of 7 per 767 trials, there is a 0.9% chance of having a click rate of 26 or more per 753 trials in the same distribution. Contextualized in the ad example, there is a .1% chance that the control ad actually performs the same as the test ad. Is the above interpretation correct? Does the test and interpretation change with the different positions (i.e. are the results of the Poisson test more appropriate for Position 3 given the small numbers)?


title: Likelihood ratio vs Bayes Factor

text: I'm rather evangelistic with regards to the use of likelihood ratios for representing the objective evidence for/against a given phenomenon. However, I recently learned that the Bayes factor serves a similar function in the context of Bayesian methods (i.e. the subjective prior is combined with the objective Bayes factor to yield an objectively updated subjective state of belief). I'm now trying to understand the computational and philosophical differences between a likelihood ratio and a Bayes factor. At the computational level, I understand that while the likelihood ratio is usually computed using the likelihoods that represent the maximum likelihood for each model's respective parameterization (either estimated by cross validation or penalized according to model complexity using AIC), apparently the Bayes factor somehow uses likelihoods that represent the likelihood of each model integrated over it's entire parameter space (i.e. not just at the MLE). How is this integration actually achieved typically? Does one really just try to calculate the likelihood at each of thousands (millions?) of random samples from the parameter space, or are there analytic methods to integrating the likelihood across the parameter space? Additionally, when computing the Bayes factor, does one apply correction for complexity (automatically via cross-validated estimation of likelihood or analytically via AIC) as one does with the likelihood ratio? Also, what are the philosophical differences between the likelihood ratio and the Bayes factor (n.b. I'm not asking about the philosophical differences between the likelihood ratio and Bayesian methods in general, but the Bayes factor as a representation of the objective evidence specifically). How would one go about characterizing the meaning of the Bayes factor as compared to the likelihood ratio?


title: Inconsistency in mixed-effects model estimation results (Stata and SPSS)

text: For my thesis there's a big chance that I will need some sort of mixed-effects specification. I have some (non-syntax) experience with SPSS but feel that it won't suffice for my analysis. I have very basic knowledge in Stata and decided to experiment more with that package. I decided to try and replicate results from SPSS in Stata for a basic model. I have data on 4059 students in 65 schools, investigating the influence of entry level score (`standlrt`) of students on their final exam score (`normexam`). In a previously followed course which had a brief introduction to multilevel modeling, my teacher provided me with a syntax in SPSS.               MIXED       normexam  WITH standlrt       /FIXED = standlrt        /PRINT = SOLUTION TESTCOV      /RANDOM INTERCEPT standlrt  | SUBJECT(school) COVTYPE(VC) .      Now I tried replicating these results in Stata but the results are not consistent. Magnitude and sometimes even sign of the betas differ. First I use `xtset school` to indicate that my data is clustered. Then I use `xtmixed normexam standlrt || school: standlrt` . What may be the cause of these inconsistent results? Thanks in advance! ps. this is not homework, and I hope I specified my first question properly. pps. a possibility may be that the 'problem' has multiple optima but I don't think this is the case in such a basic model, also because it's an uni-variate regression. Also, the iterative procedures performed while estimating may have different results, but I only think this would have big effects like sign changes. EDIT This is my Stata output               xtmixed normexam standlrt || school: standlrt          Mixed-effects REML regression                   Number of obs      =      4059     Group variable: school                          Number of groups   =        65                                                          Obs per group: min =         2                                                                    avg =      62.4                                                                    max =       198                                                               Wald chi2(1)       =    768.21     Log restricted-likelihood = -4667.8385          Prob > chi2        =    0.0000          ------------------------------------------------------------------------------         normexam |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------         standlrt |   .5570213   .0200971    27.72   0.000     .5176317    .5964108            _cons |  -.0080944   .0400842    -0.20   0.840     -.086658    .0704691     ------------------------------------------------------------------------------          ------------------------------------------------------------------------------       Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]     -----------------------------+------------------------------------------------     school: Independent          |                     sd(standlrt) |   .1214197   .0191066      .0891958    .1652852                        sd(_cons) |   .3032317   .0309434      .2482638    .3703699     -----------------------------+------------------------------------------------                     sd(Residual) |   .7440605   .0083943      .7277885    .7606962     ------------------------------------------------------------------------------     LR test vs. linear regression:       chi2(2) =   438.60   Prob > chi2 = 0.0000      And this is my SPSS output               -2 Restricted Log Likelihood    9335,677          Type III Tests of Fixed Effects(a)     |---------|------------|--------------|-------|----|     |Source   |Numerator df|Denominator df|F      |Sig.|     |---------|------------|--------------|-------|----|     |Intercept|1           |60,466        |,041   |,841|     |---------|------------|--------------|-------|----|     |standlrt |1           |56,936        |768,207|,000|     |---------|------------|--------------|-------|----|     a. Dependent Variable: normexam = final exam scores.               Estimates of Fixed Effects(a)     |---------|--------|----------|------|------|----|-----------------------------------|     |Parameter|Estimate|Std. Error|df    |t     |Sig.|95% Confidence Interval            |     |         |        |          |      |      |    |-----------------------|-----------|     |         |        |          |      |      |    |Lower Bound            |Upper Bound|     |---------|--------|----------|------|------|----|-----------------------|-----------|     |Intercept|-,008094|,040084   |60,466|-,202 |,841|-,088262               |,072073    |     |---------|--------|----------|------|------|----|-----------------------|-----------|     |standlrt |,557021 |,020097   |56,936|27,717|,000|,516777                |,597266    |     |---------|--------|----------|------|------|----|-----------------------|-----------|     a. Dependent Variable: normexam = final exam scores.               Covariance Parameters          Estimates of Covariance Parameters(a)     |-------------------------------------|--------|----------|------|----|-----------------------------------|     |Parameter                            |Estimate|Std. Error|Wald Z|Sig.|95% Confidence Interval            |     |                            |--------|        |          |      |    |-----------------------|-----------|     |                                     |        |          |      |    |Lower Bound            |Upper Bound|     |----------------------------|--------|--------|----------|------|----|-----------------------|-----------|     |Residual                             |,553626 |,012492   |44,319|,000|,529676                |,578659    |     |----------------------------|--------|--------|----------|------|----|-----------------------|-----------|     |Intercept [subject = school]|Variance|,091949 |,018766   |4,900 |,000|,061635                |,137174    |     |----------------------------|--------|--------|----------|------|----|-----------------------|-----------|     |standlrt [subject = school] |Variance|,014743 |,004640   |3,177 |,001|,007956                |,027319    |     |----------------------------|--------|--------|----------|------|----|-----------------------|-----------|     a. Dependent Variable: normexam = final exam scores.      As you can see, the log likelihoods are the same. Additionally, the fixed effects tables are the same. However the random effects are different. I'm not very skilled in interpretation yet but the results seem to differ. These are the settings for the variance-covariance matrix                Model           covariance(vartype)    variance-covariance structure of the random                                    effects              vartype                  Description             -------------------------------------------------------------------------             independent              one variance parameter per random effect, all                                        covariances zero; the default unless a factor                                        variable is specified             exchangeable             equal variances for random effects, and one                                        common pairwise covariance             identity                 equal variances for random effects, all                                        covariances zero; the default for factor                                        variables             unstructured             all variances and covariances distinctly                                        estimated      And I read online that `COVTYPE(VC) requests the default (variance component) structure for random effects, which assumes all random effects are independent.`


title: Estimated Marginal Means or means from descriptives?

text: I'm currently looking at a 2x2 Repeated Measure Anova for an undergrad project but I'm struggling to find information about which means to report. I will be reporting the Within-Subjects Anova,the Interaction and the Between Subjects Anova. I have a Profile Plot to show the Estimated Marginal Means of Measure1 which I will be using to show that there is no interaction between the two. However - I'm am struggling to find information regarding the reporting of means. I have a table of descriptive with the Mean and SD. Then towards the end of the SPSS output I have a section with Estimated Marginal Means. Which of these would I report and how would they be reported? Would this before or after details of the Anova are reported? Any help appreciated as I have been sat here for hours struggling with this.


title: Probit estimation - insignificant coefficients with a significant difference

text: I'm a postgraduate economics student but not of the statistically inclined sort :) Im running a regression to estimate the likelihood of a patent being litigated. I have two dummies for the patent owner type; let's call them `SmallFirm` and `LargeFirm`. The reference category is individuals. What I'm actually interested in is if the difference between the coefficients on Large and Small is negative, which would suggest in my particular case that small firms face a higher litigation risk. I have quite a small sample size and find each coefficient to be insignificant at the 10% level. However, when I run an test of equality of coefficients in stata (`test LargeFirm=SmallFirm`) I get that they are significantly different at the 5 % level.               Dummy     |coeff     |std err  |z-stat |P-value|       CI     SmallFirm |-.1636389 |.5219826 | -0.31 |0.754  | (-1.186706,    .8594281)     LargeFirm |-.7599469 |.533567  | -1.42 |0.154  | (-1.805719,    .2858252)          . test SmallFirm = LargeFirm          ( 1)  [Litigated]SmallFirm - [Litigated]LargeFirm = 0                 chi2(  1) =    5.88          Prob > chi2 =    0.0153      What can I conclude from this? Can I say that the significant difference does imply that small firms are associated with a higher likelihood of litigation despite the individual dummies being insignificant? Thanks so much for the help! Update! I think I got it now! Thanks Andrea: you got me on the right track. The dummies are insignificant precisely b/c my choice of reference group: neither large or small firms are statistically significantly more likely to be litigated - than individuals-. this is not what I was interested in anyway. If I run a regression with the SmallFirm dummy and an Individual dummy making LargeFirm the reference group: surely enough I do find a significant coefficient on SmallFirm! So thanks!


title: Does it make sense to consider non-binary logit?

text: Consider the nested logistic regression model’s following interpretation. A person repeatedly chooses between two different options. These options feature different levels of the same features (e.g. car characteristics). In some cases, the person will choose option 1, in some cases the person will choose option 2. The standard formula of the logistic regression $$P(Y=1|X=x_{i}) = \frac{e^{\beta x_{i}}}{1-e^{\beta x_{i}}}$$ could be interpreted such that the person evaluates options according to an utility function $U(x_{i})$, assigning utility components $\beta$ to the characteristics of the options, and choosing the option which maximizes utility. In a logistic regression we would thus estimate the implicit values of the characteristics $x_{i}$. Now suppose the setting would be such that the person does not make a discrete choice. The person could, for example, say that both options are equally attractive, that option 1 is much more attractive than option 2, or that option 2 is just slightly more attractive than option 1 etc. We would thus compute a logit specification model with the dependant variable not being 0 or 1, but within the interval $(0,1)$. Is it mathematically possible to do that? I believe it is.* If you take a look at the ML-process of calculating a logistic regression, it is nowhere required for y to be dichotomous, nor is it in fact required to be in the interval of $(0,1)$. You can fit regression models with the ML-method for any objective functions you desire, question is if it will make sense or not. It won’t exactly be logistic regression anymore, but it might be something different, similar. So the question remains: **Does that make any sense?** Is it valid to think of a binary choice problem the way I portrayed it, and to make use of additional information on relative weighting of preferences (if available) by assuming y \e (0,1)? Would that make better use of available information, as in the case of y ~ 0.5 we know that the valuation of the options is similar?


title: Why do I get 100% error rate in unsupervised random forest, and how does the unsupervised pattern work in "randomForest" R package

text: I tried to use random forest to classify microarray data. Basing on research of L.Breiman and Tao Shi, I constructed a synthetic data base using bootstrap methods (Assuming it is a matrix with samples on row and genes on column, for each gene in each sample, values are selected with replacement in the gene column, descripted in http://www.stat.berkeley.edu/~breiman/RandomForests/cc_home.htm#unsup ). The original data is labeled as 1, and the synthetic data is labeled as 2. The combined data is used as input for random forest as supervised data. The wired stuff is that I got a error rate of 100%. In my view, the worst case would be the error is 50%, which means randomly participate the combined dataset. I tried this approach for many times, it always gave nearly 100% error rate. In practical, my data set is 7*297 matrix, 7 samples and 297 genes. Besides, I have a problem in understanding how unsupervised pattern work in the R "randomForest" package, and how it calculate the proximity. From the source code, I get this:               if (!is.null(y)) {         if (length(y) != n)              stop("length of response must be the same as predictors")         addclass <- FALSE     }     else {         if (!addclass)              addclass <- TRUE         y <- factor(c(rep(1, n), rep(2, n)))         x <- rbind(x, x)     }      It seems that the new data set is treated as combined raw data. However, I didn't get observed result if the assumption is true:               iris.urf <- randomForest(iris[, -5])     iris.urf$proximity[1:5,1:5]     MDSplot(iris.urf, iris$Species)          xx.1 <- randomForest(as.matrix(rbind(iris[, -5],iris[, -5])),factor(rep(1:2,each=150)),proximity=TRUE)     xx.1$proximity[1:5,1:5]     xx.2 <- xx.1$proximity     xx.1$proximity <- xx.2[1:150,1:150]     MDSplot(xx.1, iris$Species)      The proximity is very different, but the mdsplot gives similar pattern. Could someone help me on these?


title: Origin of Z-estimation

text: Van der Vaart [1] uses the term Z-estimation for the estimating equations approach. Is this the first use of the term? [1] Vaart, A. W. Van Der. Asymptotic Statistics. Cambridge University Press, 2000.


title: What is the difference among Indicator, Index, Variable and Measure?

text: What is exactly the difference among **Indicator, Index, Variable and Measure**? I would appreciate some reference.


title: What time-series model to use for mostly linear, but limited values?

text: I'm interested in fitting a model to several types of time series that all behave similarly:   * Depth of a work queue   * Disk space available   * Idle CPU time available None of these values can go below zero, but otherwise their behavior is relatively independent of value. (e.g. disk space is as likely to move from 25% to 26% as it is from 75% to 76%.) Some of these, like disk space and idle CPU time, also have upper limits that they cannot pass - a disk cannot be more than 100% full. As long as these values stay away from their limits, I suspect traditional time series models like ARIMA or Holt-Winters would model them fairly well, but I'd like to encode the behavior at the limits. What models can encode these behaviors?


title: NMDS for biomass

text: I would like to make a NMDS with biomass of different prey groups in stomach content of fish. I have already made one where the data matrix consists of 0 and 1, and this one went fine but are not able to do it when adding the actual weights. I have used the MASS and vegan package. When using binomial data matrix:               mymetaMDS <- metaMDS(comm, distance = "jaccard", na.rm = TRUE, binomial = TRUE)       and this works fine. But when I try to create the dissimilarity matrix with the weight data matrix:               mymetaMDS <- metaMDS(comm, distance = "jaccard", na.rm = TRUE, binomial = FALSE)       > Error in if (any(dist < -sqrt(.Machine$double.eps))) warning("some > dissimilarities are negative -- is this intentional?") : missing value where > TRUE/FALSE needed In addition: Warning messages: 1: In distfun(comm, method > = distance, ...) : you have empty rows: their dissimilarities may be > meaningless in method "jaccard" 2: In distfun(comm, method = distance, ...) > : missing values in results Any suggestions any one?


title: Are all Levy processes memoryless?

text: We know that the two canonical Levy processes, namely the Wiener process and Poisson process, are both memoryless, so I wonder if there are any Levy process that is not memoryless. Specifically, are these two properties equivalent: independent stationary increments <=> memorylessness ? Thanks.


title: Does non-zero correlation imply dependence?

text: We know of the fact that zero correlation does not imply independence. I am interested in whether a non-zero correlation implies dependence - i.e. if $\text{Corr}(X,Y)\ne0$ for some random variables $X$ and $Y$, can we say in general that $f_{X,Y}(x,y) \ne f_X(x) f_Y(y)$?


title: Can a regression coefficient be explained both ways in logistic regression?

text: In logistic regression, depending on the sign of the coefficient we can say, for example, that when predicting political interest, a coefficient of $1.226$ for following the news indicates that people are more likely to be interested in politics when they are following the news. However, can we say that based on this coefficient, the more people are political interested the more likely they are to follow the news? In other words, does it work both ways?


title: How I can deal with too many variables in training a data set?

text: I am trying to train a predictive model on whether a given person is ( male or female) based on behavior cues we've obtained from online surveys.   1. The dependant variable will be a binary ( 1 or 0 ).    2. For the training data set, I have a list of 150+ numeric variables ranging from issues as far apart as "political leanings" to "amount daily exercise". These are all indexed numbers and range from 1 to 100.  So my problem is, I have too many variables! I don't really want to assume some variables are more important than others because there may be something I don't know about. And I don't want to throw away important correlations between variables. But how can I build a good, strong, predictive model using all the data I have?


title: How do you prove the correlation of two linearly related variables is equal to one

text: Assume that $X$ is a zero mean random variable and $a$ and $b$ are constants. Show that if $Y = aX + b$ for some $a > 0$ than correlation $(X,Y) = 1$.


title: Sample size for one sample repeated measures

text: I am conducting a study in which I want to follow a group of subjects over time. My outcome of interest is continuous and I will measure it at 4 time points. I am trying to calculate the sample size necessary but all The sample size formulae i am finding involve two or more groups. The only table I found for a one group repeated measures design is by bracikowski. Can I rely on it for my sample size calculation? Any other recommendations? Thank you


title: Measure the relation/association between the outcome and the independent variables

text: I am running a model (logistic regression) with 20 independent variables in R. Before running the model I calculated the correlation between all the variables and finally selected my variables by also checking "visually" the histograms of each variable in the case of presence and again in the case of absence. In situations where I don't see any obvious distribution associated to both presence & absence, I discard the variable. I would like to make "official" calculations for the level of relation between Presence/Absence and each variable (how much each variable contributes to the Presence/Absence), for example with `Cramer's V index`, but the available function I find is from the package `vcd` and has some limitations: doesn't give the `Cramer's V` (as well as the Phi-Coefficient Contingency Coeff.) for each independent variable, and it doesn't run for one independent variable. I might be missing some other obvious way to do this. Any help is appreciated.


title: Logistic/Association model - numerical vs. categorical covariates & target

text: I am using a powerful analytics tool built by Oracle (Data Miner). I am trying to determine the probability of death for a patient given the presence of certain diagnosis codes. There will be up to 40 possible codes as column references (covariates) in the model. Any reason why I must convert the presence indicator to a numeric 0/1 design variable ? Same thing with the target variable (death) - must it be 0/1 (0=alive, 1=death) ? I noticed that SAS users do this regularly, but my modelling package allows "Y","N" instead. Also, I can have a value of Null which indicates an unknown.


title: Is two years enough for panel data analysis?

text: I have around 800 companies for only two years period. However, around 200 of them have only one year observation. Is it still possible to conduct panel data analysis with such data Thank you


title: How to define the termination condition for gradient descent?

text: Actually, I wanted to ask you how can I define the terminating condition for gradient descent. Can I stop it based upon the number of iterations, i.e. considering parameter values for, say, 100 iterations? Or should I wait such that the different in the two parameters values 'new' and 'old' is very small to the order of lets say $10^{-6}$? This will definitely take much time. What is the best way? In my case even one iteration takes significant time. In this situation if I wait for the 2nd condition it might even take weeks I guess. So which approach should I use. How to tackle this scenario?


title: Differences between poisson.test and E-test when testing Poisson parameters

text: Suppose we have two independent Poisson-distributed variables $X_1$ and $X_2$. We want to test whether the Poisson parameters are equal, i.e. whether $\lambda_1=\lambda_2$. Now we have 4 distinct statistical _exact_ tests to choose:   1. E-test (see Krishnamoorthy and Thomson, A more powerful test for comparing two Poisson means, _Journal of Statistical Planning and Inference_ 119 (2004) 23–35; see also Checking if two Poisson samples have the same mean on Cross Validated)   2. `poisson.exact(tsmethod="central")`   3. `poisson.exact(tsmethod="minlike")`   4. `poisson.exact(tsmethod="balker")` (from exactci R package) Now, given all those tests are labeled as "exact", one would expect all yield the same p-values. Contrary to this, the quoted paper clearly illustrates that methods 2-4 give different significance. Furthermore, I personally implemented the E-test and found that this test gives yet another, distinct result. Why is that?


title: What is the conventional definition of recurrence-free survival?

text: Hello to all the biostatisticians and epidemiologists out there, I have looked up and down for a standard definition of recurrence-free survival, and the issue I'm having is determining if the standard includes DEATH as an event or not. I have seen some studies censor at death and others include death along with recurrences as events. If death is traditionally included as an event, what is the rationale for that? It certainly is not a recurrence. Also, is there a difference between relapse-free survival and recurrence-free survival? Thanks for your help!


title: SAS ARIMA: estimate without identify? or auto.arima in SAS?

text: I have several groups of data for time series analysis, and I need some automatic way to do the analysis. I try to find a way in SAS to automatically do the job, but I can't find one corresponding to `auto.arima()` in R. I have two questions regarding this:   1. The closest one is to use _identify_ with _scan_ and _esacf_ to find tentative orders for arima model. However, that doesn't test unitroot or whether the data is just a whitenoise. Also, I still need to manually figure out the orders. **Does anyone know if there is any procedure in SAS that can do the analysis automatically?**   2. I tried to write a macro to do the job automatically. I need to fist run _identify_ to check unitroot, whitenoise, orders etc. and based on the results, I need to run _estimate_ accordingly with corresponding orders. However, when I write two separate _proc arima_ procedures, SAS requires me to run _identify_ again in the second procedure, otherwise the _estimate_ can't be executed. For example:              PROC ARIMA data = xxx;          IDENTIFY VAR = yy scan esacf;     RUN;     /* other code */          PROC ARIMA data = xxx;          ESTIMATE p=1 q=1; <---------------- ERROR: must run IDENTIFY first     RUN;      **Is there a way to use the results of the _identify_ in the first ARIMA for the _estimate_ in the second ARIMA procedure?** Otherwise, it's just a waste of time to run the _identify_ twice.


title: How to carry out repeated time-to-event analysis in SPSS?

text: **Please note that I also posted this same question in an SPSS specific forum, however if anyone here can also help, I would appreciate it.** I have a question regarding repeating time-to-events measurements and after a quick search I really did not find anything that directly answered my question. So here goes... Context (an example of the model I am trying to create): I have say 25 patients who have all undergone the same procedure and who meet eligibility criteria for a study. Date of discharge from hospital is "Day 0" and I am tracking multiple events following their discharge that cause these patients to be readmitted. The type of event is not necessarily important at this stage (could be anything from infection, to bleeding, to device failure, etc.), as any event causing readmission is classified as `1=true` or `0=false`. Patients who then receive a subsequent procedure have reached their end-point are censored from that point on. The issue I am having is that each patient may have multiple readmissions following treatment, thus:               Patient 1: day 0-->6* (readmission), day 7-->25* (readmission), day 26 -->65* (readmission);      Patient 2: day 0-->34* (readmission), day 35-->38* (readmission)     Patient 3: day 0-->50*(readmission), day 51-->53**(subsequent procedure), *Day 54-->60(readmission)*     ..etc.       where `*` notes the day being readmitted to hospital. and text following the 'double asterisks' means any readmissions at this point forth are excluded because they have reached end-point. Incidence per patient-year at say day 30 was just a matter of calculation by hand, **but if** I wanted to perform some sort of Kaplan-Meier survival type estimate or a Cox Hazard regression for risk of readmission using $x$ and $y$ covariates, is that possible to do in SPSS? Even though with 25 patients, the study may be underpowered to come to any meaningful conclusion with further analysis beyond descriptives, is this possible and how? Thus far, my previous works have been limited to Kaplan-Meier survival, log- rank analysis, and Cox-regressions with single definite event (i.e., death) that is continuous and non-repeating... so I'm treading into new water. If anyone has any advice or can point to me to specific reading materials, that would be greatly appreciated.


title: Help in viewing the output structure of a neural network

text: 0 down vote favorite 1 I used nftool in Matlab 2012 and trained a network. I gave the training inputs as x=[250:1] and targets as t=[250:1]. I used 10 hidden layers. I trained the network and got the results and exported the data to workspace. I finally got net as the variable in the workspace. I want to realize this network in hardware. So, I want the exact weights and bias values. I searched the net and typed net.IW and got 10 values, net.LW and got yet another 10 values and net.b and got 11 values. My network diagram is as shown below: Output Neural Net ![enter image description here](http://i.stack.imgur.com/hx1kG.png) If IW is the Input Weight and LW is the Layer Weight, I realized the network diagram as shown below: My realization of output Neural Network ![enter image description here](http://i.stack.imgur.com/xNywy.png) I got 10 values for IW, which would fit into 10 hidden layers, and 10 values for LW, which will fit into 10 hidden layers according to my realization. I got 11 values for bias, which will fit in to my network, since there are 11 (b) blocks in the network. But, i am missing one IW and one LW value. I want to know whether there is a mistake in my realization of the output of the network or I missed any of the values. Please help.


title: How is the kurtosis of a distribution related to the geometry of the density function?

text: The kurtosis is to measure the peakedness and flatness of a distribution. The density function of the distribution, if it exists, can be viewed as a curve, and has geometric features (such as curvature, convexity, ...) related to its shape. So I wonder whether the kurtosis of a distribution is related to some geometric features of the density function, which can explain the geometric meaning of kurtosis?


title: What is the intuitive sense behind the purpose and mechanics of Sufficient Statistics?

text: The definition of a sufficient statistic is: Let $X_1,...,X_n$ be a random sample from a distribution indexed by a parameter $\theta$. Let $T$ be a statistic. Suppose that, for every $\theta$ and every possible value $t$ of $T$, the conditional joint distribution of $X_1,...,X_n$ given that $T=t$ depends only on $t$ but not on $\theta$. Then, $T$ is a sufficient statistic for parameter $\theta$. I feel like I know several pieces of the puzzle (like the factorization theorem) to understanding sufficient statistics but do not have the overall theory down. My main questions are: 1) Why do they say that $T$ is a sufficient statistic **_for parameter $\theta$_**? If $\theta$ were the population mean of a normal distribution, say $\mu$, does it mean that anytime we want to find the probability of, say, $X_1,...,X_n$ occurring in a certain way, that we don't need the value of the mean of the population? 2) In real-life, why do we want to use a sufficient statistic? Is seems that just calculating the statistic shouldn't be that much work (such as the sum of X's) so why do we need it? Thanks!


title: Log likelihood - understand depper

text: I want to use log likelihood formula to relate between two items. The formula is:               LLR = 2 sum(k) (H(k) - H(rowSums(k)) - H(colSums(k)))      When this is the table: > ...............................| **Event A............................** | > **........ Everything but A** > > **Event B.................** | A and B together (k_11) .| B, but not A > (k_12) > > **Everything but B** | A without B (k_21)...........| Neither A nor B (k_22) where H is Shannon's entropy, computed as the sum of (k_ij / sum(k)) log (k_ij / sum(k)). I know it’s trying to assess how unlikely it is that the overlap between the two items is just due to chance. But I don't understand the logic - why `H(k) - H(rowSums(k)) - H(colSums(k)` will lead this answer? what is the logic behind this?


title: Feature-based malware detection

text: My aim is to detect a malware by testing it with the extracted features. I have extracted the binary features of 150 malware and 150 benign files by using the hex editor tool for unpacking the files and net beans for generating a dataset into an Excel sheet. Now I have features of malware and benign datasets in separate sheets. Now I want to test a new file with the existing features to say whether it is a malware file or a benign file, for this I'm working with WEKA tool. can I achieve my goal with this tool? If not, Which tool is better for my project. I just want to compare the new file automatically with those two malware and benign datasets when given as input and should say if it is malware file or benign file. Thank you...


title: Question about Chow forecast test

text: I would like to check my understanding of the Chow Forecast test.The test is used to perform a test on the following models: $y_1=\beta X+\epsilon_1$ and $y_2=\beta X+ \nu + \epsilon_2$, with $H_0:\nu=0$ and $H_1: \nu\neq 0$. Now I wonder whether that means that the following holds true: \begin{align*} F&=\frac{(e_R'e_R-e_1'e_1)/g}{e_1'e_1/(n-k)}\\ &=\frac{\left(\sum\limits_{i=1}^{n+g} (y_i-x_i'\beta)^2-\sum\limits_{i=1}^n (y_i-x_i'\beta)^2 \right)/g}{\sum\limits_{i=1}^n (y_i-x_i'\beta)/(n-k)}\\ &=\frac{\sum\limits_{i=n}^{n+g} (y_i-x_i'\beta)^2/g}{\sum\limits_{i=1}^n (y_i- x_i'\beta)/(n-k)} \end{align*} Are these true, and if not, why not?


title: Fitting of bivariate data to a self-defined probability density function

text: I have a bivariate set of data points which I want to fit to a self-defined distribution (i.e. not standard normal or chi-square or like that, a different, let's say "new" density function). I would wish to fit it via method of maximum likelihood, but my density function has 2 infinite sums in it so strictly speaking I don't have a closed form expression. The density function is also quite a clumsy expression. I know about the "nlm" function in R, but I have a very limited knowledge about R. I am very familiar with SAS and fairly familiar with Mathematica. Does anyone have any suggestions as to how I can proceed? Thus, for a given set of bivariate data and a self-defined bivariate density function, I wish to estimate 5 parameters via the method of maximum likelihood. Hoping someone can help me :)


title: Convert z-scores to percentile rank

text: Does anyone know the formula to convert a z-score to a percentile rank instead of using the z-score table? I searched the website, but did not find the formula. thank you


title: How to judge a nonparametric estimator empirically?

text: Suppose we have an estimator $\hat{f}(Y_1, \cdots, Y_n)(t)$ for the fixed design nonparametric regression model: $$ Y_i = f(\frac{i}{n}) + \epsilon_i, \; f \in L^p[0,1]. $$ The estimator at hand might have very good theoretical properties (e.g. achieves certain oracular risk up to a multiplicative factor/is asymptotically minimax/etc.). But how does one translates this empirically? I am guessing something along the following lines: fix $f \in L^2$ and let's say we use the $L^2$-loss. Compute from simulation the quantity $$ \frac{\|\hat{f} - f\|^2_{L^2}}{ \|f\|^2_{L^2}}, $$ so basically the MISE divided by the energy of the signal. If this ratio is "small" in some sense...we're happy? It looks like the signal-to-noise ratio ($sd$ denotes empirical standard deviation) $$ \frac{sd(f(\frac{i}{n}))}{sd(Y_i)} $$ should also be taken into account somehow.


title: accessing learned coefficeints (coef_) of logistic regression in sklearn

text: I am wondering about how to interpret the coef_ variable in the logistic regression class of sklearn. Given a dataset with m features and n categories, coef_ seems to be a matrix with the size of (n-1, m) -- In Python, it is a list (size n-1) of a list (size m). In 2 category classification (i.e., the used categories are 0 or 1), coef_ definitely seems to be the coefficients for computing the likelihood of category 1. In multi-category classification (where the number of the categories is more than 3), how does each list of the weights in coef_ correspond to a category? For example, if 0, 1, and 2 are used as a category, does the first list in coef_ correspond to the category 0 or to the category 1?


title: Unsure how to calculate mean square error of a variable with a joint distribution

text: I'm working through a homework assignment and I'm having an inordinately difficult time trying to calculate the MSE of two different estimators. The homework is described as this: Let $X$ and $Y$ be random variables with joint density: $p(x,y) = 2$ if $x\in A$ It is $0$ otherwise. $A$ is a set surrounded by $X = 0, Y = 0$, and $X + Y = 1$. I have to calculate the MSE of 2 estimators: \begin{equation} \hat{X} = E[X|Y] \end{equation} and \begin{equation} \hat{X} = Y \end{equation} I've seen the proof that $E[X|Y]$ is the better estimator, but I am absolutely stumped on how to actually calculate the MSE $E[(\hat X - X)^2]$ for either estimator. I've spent hours reading and trying to figure this out but I have almost no experience with statistics and it is killing me in this problem. Any help would be appreciated. Edit: $E[X|Y] = \int_0^1 \int_0^{1-y} 2xy dx dy = 1/12$ $E[Y] = \int_0^1 2ydy = 1$


title: Why do we estimate the mean response in Confidence interval but predict individual outcome in prediction?

text: I have tried to understand the differences. Just to be clear, I _think_ I understand the two sources of variation in prediction comes from the variation in the distribution of the _location_ of Y and _within_ variation of Y. The thing I do not get is, _why_ we do not use the estimate of the mean response (like in estimation) but use the individual outcome of Y (thus the two variation) in prediction? Or to put it other way round, why in estimation we did not consider variation in the distribution of the location of Y? Is it because, Y is a rv in prediction? If so, why Y is not rv in estimation?


title: How to use Kernel Density Estimation for Prediction?

text: I would like to apply KDE to inventory replenishment, but I am not sure how to use the analysis to predict future sales based on past sales. Given a set of data and having applied KDE to it (probably using a Gaussian distribution), how do I make a prediction about the future? Thanks for any help! Please let me know if I can clarify - I'm only starting to pick up the language for talking about KDE ... I'm glad to do reading on my own - pointers to any resources would be welcome.


title: Conditional expectation subscript notation

text: This should be a relatively simple question. I'm trying to confirm my understanding of the subscript notation on expectations when the subscript denotes a conditioning. In the example $$E_{Y|X}[(Y-f(X))^2|X]$$ the subscript denotes the distribution over which you take the expectation, so we would want to use $p_{Y|X}(y|X)=P(Y=y|X)$ as the distribution in our expectation and sum over the density of $y$. And the argument $$[(Y-f(X))^2|X]$$ means that when we take the expectation we should consider the value of $X$ in $(Y-f(X))^2$ to be given. This would imply that the statement can be rewritten: $$ E_{Y|X}[(Y-f(X))^2|X]=\int_{-\infty}^{\infty} [y-f(X)]^2 p_{Y|X}(y|X) dy. $$ And if I wanted to find $E_{X}E_{Y|X}[(Y-f(X))^2|X]$ I could rewrite that as: $$ E_{X}E_{Y|X}[(Y-f(X))^2|X]=\int_{-\infty}^{\infty} \left( \int_{-\infty}^{\infty} [y-f(X)]^2 p_{Y|X}(y|X) dy]) \right) p_{X}(x)dx. $$ Is this all correct? **Correction** Updated last line to fix an error that was caught in the comments: $$ E_{X}E_{Y|X}[(Y-f(X))^2|X]=\int_{-\infty}^{\infty} \left( \int_{-\infty}^{\infty} [y-f(x)]^2 p_{Y|X}(y|x) dy]) \right) p_{X}(x)dx. $$


title: Multinomial regression analyses

text: My dependent variable is categorical (with 3 levels) and my predictor variables are a mix of continuous (age in months, test 1 score, test 2 score, test 3 score) and categorical (gender). I believe I should run a multinomial regression, but when I do the results are really uninterpretable because each age is basically treated as a category (and the same for the test scores). Does anyone know a different test that would be appropriate to use in my case? I would really appreciate any help.


title: Use variables selected from DISTLM multivariate regressions in CAP

text: I would like to know if the following analysis strategy is correct conceptually, if it makes sense. I found difference in fish assemblage structure between two estuarine sectors, through PERMANOVA. I would like to identify what environmental factors are associated with this difference. Then, I used DISTLM (Distance‐based linear modelling) utilising information criteria to find the best model of variables related to the multivariate data cloud. I then run a canonical analysis of principal coordinates (CAP) to find axes maximizing this spatial differences with linear combinations of variables from the best DISTLM model. I'm not so sure if I can use these variables from DISTLM in the CAP. If not, can anyone suggest an alternative analysis strategy?


title: Derivation of the equation of ridge regression

text: I am having some issues with the derivation of the equation of ridge regression I know the part without the regularization param like $\beta = (X^TX)^{-1}X^Ty$ But after adding the L2 term $\lambda||\beta||_2^2$ how come the solution is $\beta = (X^TX + \lambda I)^{-1}X^Ty$ I mean how come the gradient of L2 term is $\lambda I$


title: Zero to one proportion data, with quadratic term?

text: I have a data set with the following: N = 60; x = developmental stage (range 25 to 44); y = proportion of 10 minute trial performing a behavior (range 0 to 0.81; 30 zeros) A scatterplot produces a quadratic looking curve where those in mid- development clearly performed the behavior for more time. Most of the zeros are in the youngest and oldest individuals. If I break up the data into 5 groups according to developmental stage, an ANOVA/Tukey strongly supports this pattern. However, I would like to analyze this data continuously without breaking it into groups. I have considered arcsine square root transformed proportion data in a linear regression, but I am unsure if that can incorporate a quadratic term, and this analysis results in a very small R squared value (less than 0.1). I have also considered arcsine square root transformed proportion data in a GLM containing a quadratic term or a beta regression (zeros??), but am not sure where to go from here. I am planning to say in the paper that the individuals in mid-development perform the behavior more than those in early or late development, but am struggling to interpret the data in a way that supports that statement. I appreciate any suggestions, thank you!


title: Creating formula object for coxph()

text: I am trying to create the formula object for `coxph()` manually, as I want to reproduce it in RSRuby. Note that the response variable needs to be a `Surv` object, the output of the `Surv()` method. The `Surv` object is basically a matrix with columns `time` and `status`. Now, I am doing the following:               length_of_stay <- c(13, 4, 5, 21, 33, 10)      exited_care <- c(1, 0, 1, 1, 0, 1)      survival_object <- Surv(length_of_stay, exited_care)     gender <- c(1, 0, 1, 0, 1, 0)     mydata <- data.frame(length_of_stay = length_of_stay,                           exited_care = exited_care,                           gender = gender)     data_frame_for_formula <- data.frame(cbind(survival_object, gender))     my_formula <- formula(data_frame_for_formula)     mod.los <- coxph(formula = my_formula, data = mydata)      and getting the following error message:               "Error in eval(expr, envir, enclos) : object 'status' not found"      `traceback()` shows the formula is               time ~ status + gender      How can I avoid that and make the survival object the response variable?


title: How to decide whether to reuse old code or reinvent the wheel?

text: For a long time now, I have been thinking about working with neural networks and genetic algorithms. I have never been able to decide whether it makes sense to start writing my own code, or to reuse the many, many options made by others. In my experience, assuming the programmer knows what he or she is doing to a reasonable extent, there is very little that replaces C or C++ code crafted to the individual users preferences. For one thing, it empowers you to do almost anything with the code due to having complete knowledge of its workings. One can for instance incorporate whatever new advances that one happens to read about in the literature. On the other hand, one person will often find it hard to do as much as a team of people. As an individual, it's always possible that one will miss various bugs. There are many other disadvantages to being a loner. So, for those who have some experience in machine learning, I would like to hear what actually works in the real world. Is it more common to use a few well-known programming packages or do many people write their own custom implementations? I am specifically interested in the areas of RBM, neural networks, genetic programming and clustering. Assume that I am willing to put a lot of time into learning about machine learning and that I will put at least as much time and effort into it as a graduate student would.


title: Interpreting "elliptic" shape in visualization of high-dimensional dataset, e.g. COIL-20?

text: There exist an image dataset referred to as COIL-20. There are several paper that describes different dimensionality reduction methods and apply them to produce a 2d plot of this dataset. These pictures often contain one or several "elliptic" shapes. Is there any particular explanation for this? If yes, can we extend it to other datasets? E.g. if we observe the similar"elliptic" shapes on completely different dataset (not images), is it possible to infer something about its structure? For example, the picture by t-SNE: ![enter image description here](http://i.stack.imgur.com/pdbVW.png) or from here: ![enter image description here](http://i.stack.imgur.com/97iAL.png) and here: ![enter image description here](http://i.stack.imgur.com/Rw1Od.png)


title: Consider the sum of $n$ uniform distributions on $[0,1]$, or $Z_n$. Why does the cusp in the PDF of $Z_n$ disappear for $n \geq 3$?

text: I've been wondering about this one for a while; I find it a little weird how abruptly it happens. Basically, why do we need just three uniforms for $Z_n$ to smooth out like it does? And why does the smoothing-out happen so relatively quickly? $Z_2$: ![2](http://i.imgur.com/ER5qI.gif) $Z_3$: ![3](http://i.imgur.com/5y620.gif) (images shamelessly stolen from John D. Cook's blog: http://www.johndcook.com/blog/2009/02/12/sums-of-uniform-random-values/) Why doesn't it take, say, four uniforms? Or five? Or...?


title: Correlation coefficient for Boolean data

text: I have a data set consisting of a number of variables: all are Boolean true/false type. Can I simply count the R correlation coefficient between two variables substituting 0 and 1 to get correlation matrix?


title: Is there an R package for continuous time longitudinal binary responses?

text: The `bild` package appears to be an excellent package for serial binary responses. But it is for discrete time. I would like to specify a smooth function of time for the odds ratio connection of the current response Y with binary responses measured at earlier times, or at least a first-order Markov version of this. I believe this is called alternating logistic regression. Does anyone know of an R package that handles continuous time, i.e., measurement times can be at any follow-up time? I don't need random effects in the model.


title: Solve discriminant functions to get decision boundary in R

text: I am trying to solve the discriminant functions in R for two, 2-D gaussian distributions and trying to find the decision boundary. The discriminant function is given by, $g(x) = -.5x'\sigma^{-1}x + \sigma^{-1}\mu x$ where, $x'$ is $x$ transpose, $\mu$ is the mean and $\sigma$ is the covariance of the distribution. $\sigma^{-1}$ is inverse of the covariance matrix. I am trying to solve for $x$ by using, $g_1(x) - g_2(x) = 0$, where $g_1$ and $g_2$ represent the two distributions. The mean and the covariance matrices for the distributions are in R code below:               library(MASS)     mu1 = c(3,6)     mu2 = c(3, -2)          sigma1 = matrix(c(.5,0,0,2),2,2)     sigma2 = matrix(c(2,0,0,2),2,2)     c1 = mvrnorm(n, mu1, sigma1)     c2 = mvrnorm(n, mu2, sigma2)          plot(c1, xlim=range( c(-10,10) ) , ylim=range( c(-10,10) ) , col="red")     points(c2, col = "blue")      How to solve this in R to get the decision boundary between the distributions?


title: Inconsistency of findings between Kappa and ICC for IRR study

text: I used the following R script to calculate the Kappa and ICC for a measured done by two raters on 23 subjects:               library(irr)     temp <- structure(list(value.x = c(10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 2, 10, 10, 2),                             value.y = c(8, 8, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 8, 10, 10, 10, 10, 5, 10, 8, 5)),                       .Names = c("value.x", "value.y"), row.names = c(NA, -23L), class = "data.frame")          kappa2(temp[,c("value.x","value.y")],weight="squared")           Cohen's Kappa for 2 Raters (Weights: squared)           Subjects = 23         Raters = 2          Kappa = 0.121                   z = 1.61        p-value = 0.106           icc(temp[,c("value.x","value.y")],model="twoway",type="agreement",unit="average")           Average Score Intraclass Correlation             Model: twoway         Type : agreement              Subjects = 23           Raters = 2         ICC(A,2) = 0.892           F-Test, H0: r0 = 0 ; H1: r0 > 0       F(22,22.2) = 8.98 , p = 1.27e-06            95%-Confidence Interval for ICC Population Values:       0.746 < ICC < 0.954      The measurement is actually a 5-point scale which variables are `c(0,2,5,8,10)`. I am wondering why Kappa and ICC gives so much different result: Kappa concluded that there is poor agreement, but ICC concluded that there is excellent agreement. Which one should I believe?


title: Bayesian meta analysis of residual standard deviation using BUGS

text: I have a small number of studies where the same model was calculated and want to infer the typical standard residual deviation. I have the degrees of freedom, the sum of the squares SSR and the mean of the square of the residuals MSR, but not the residuals themselves. My current plan is to do a bayesian meta analysis using either WinBugs or jags using two steps: 1) Assume the residual standard deviation is really the same over all studies and just try to estimate it more precisely. Let's call the overall residual variance res.var With this, I already run into practical problems with the BUGS modelling language. I wanted to model the residual sum of squares using a chi square distribution, but have problems scaling the data accordingly. My problem is what I wanted to use is that for each study $$ \frac{SSR}{\text{res.var}} \sim \chi^{2}_{k}. $$ But I have problem putting this into a bugs model. I tried to use               model     {     for (i in 1:N)     {     chisquare[i] <- sumsq[i]/res.var     chisquare[i] ~ dchisqr(df[i])     }     res.var~dunif(0,1)     }      but this gives an multiple definitions of chisquare error. I tried transforming the first definition of chisquare so that sumsq is on the left side, but data can not be put on the left side of a logical node. Any practical advice? 2) Should I get this model to run, I would also like to allow for some variation of the residual standard deviation over studies, due to there being at least some small differences in the study population not included in the meta analysis model but I am not certain what a good hyperdistribution for the variance could look like. Any pointers?


title: Intuition on upper bound of the number of mistakes of the perceptron algorithm and how to classify different data sets as "easier" or "harder"

text: I was trying to understand the following result more intuitively (for linearly separable data): $$ k \leq \frac{R^2}{\gamma_g^2}$$ where: k = is the number of mistakes the perceptron algorithm does R = is the distance at which the data points are confined to. $\gamma_g = \text{the smallest geometric margin for the current data set i.e.} = min_{1 \leq i \leq n}\{ \gamma(i) \}$ $\gamma(i) = \text{the perpendicular geometric distance from the decision boundary for data i .}$ The way I wanted to understand better is, I was thinking of assessing the level of difficulty of a data set according to the margin they could have or the distance from the origin they had. i.e. by the number of mistakes the algorithm could make. i.e. by looking at how big/small each upper bound is different for each training data set. So for example, I currently think of this result as the following: if we can achieve a large geometric margin (i.e. if there exists a large geometric margin) then we can have higher confidence of our classifier, thus, for large margins, the denominator for the above result will be larger, which means that the number of mistakes (k) will decrease, which in turn means that, if a higher margin is possible, then the classification task is "easier". Which makes sense, if there exists a decision boundary that can clearly distinguish +'s from -'s then, the classification task should be considered "easier". However, I was having difficulties coming up with a similar explanation for the numerator. Obviously, the larger the numerator is, from the above equation, the classification task might have more mistakes (so one can consider the task "harder" for the current data set). However, I was not really sure why R has that effect on the number of mistakes nor did I have a good intuition on why R affect the perceptron the way it does. Thanks in advance! :) NOTE: I know the number of mistakes perceptron makes depends on the initialization when it starts cycling. However, that is not what I care about. I care about just arguing the intuition for the number of mistakes soley on the upper bound of the perceptron. How to initilize it is a disscussion for a different quesiton.


title: RandomForest - MDS plot interpretation

text: I used randomForest to classify 6 animal behaviours (eg. Standing, Walking, Swimming etc.) based on 8 variables (different body postures and movement). The MDSplot in the randomForest package gives me this output and have problems interpreting the result. I did a PCA on the same data and got a nice seperation between all the classes in PC1 and PC2 already, but here Dim1 and Dim2 seem to just seperate 3 behaviours. Does this mean that these three behaviours are the more dissimilar than all other behaviours (so MDS tries to find the greatest dissimilarity between variables, but not necessarily all variables in the first step)? What does the positioning of the three clusters (as e.g in Dim1 and Dim2) indicate? Since I'm rather new to R I also have problems plotting a legend to this plot (however I have an idea what the different colours mean), but maybe somebody could help? Thanks a lot!! ![RandomForest MDSplot of 6 different animal behaviours](http://i.stack.imgur.com/xiwVU.jpg) I add a plot made with the ClassCenter function in RandomForest. This function also uses the proximity matrix (same as in the MDS Plot) for plotting the prototypes. But just from looking at the datapoints for the six different behaviours, I can't understand why the proximity matrix would plot my prototypes as it does. I also tried the classcenter function with the iris data and it works. But it seems like it doesn't work for my data... Here is the code I used for this plot be.rf <- randomForest(Behaviour~., data=be, prox=TRUE, importance=TRUE) class1 <- classCenter(be[,-1], be[,1], be.rf$prox) Protoplot <- plot(be[,4], be[,7], pch=21, xlab=names(be)[4], ylab=names(be)[7], bg=c("red", "green", "blue", "yellow", "turquoise", "orange") [as.numeric(factor(be$Behaviour))]) points(class1[,4], class1[,7], pch=21, cex=2, bg=c("red", "green", "blue", "yellow", "turquoise", "orange")) My class column is the first one, followed by 8 predictors. I plotted two of the best predictor variables as x and y. ![RandomForest ClassCenter plot for 6 different animal behaviours](http://i.stack.imgur.com/lYsZe.jpg)


title: Report coefficients of a factor in a GLM

text: I have a GLM and my explanatory variables are a factor of 9 levels (X1) and a numeric variable (X2). If I have to report the formula: $Y = B_0 + B_1X_1 + B_2X_2$ I have 1 coefficient for the intercept and 1 coefficient for the numeric variable, that's ok, but I have 8 coefficients for the factor, how do I report them? I mean, Bo is one number, B2 is one number, but B1? do I sum the eight numbers? or I report the mean of the 8 numbers, or I have to report the formula with the eight numbers? Thank you.


title: Density calculation in K Means clustering

text: With the K Clusters generated using K Means Clustering, how do we calculate the density of each cluster? Is there any formula for it?


title: Forecasting a complex time series by splitting into subseries

text: I have finance data that I need to forecast out for 7 years. My data is generally debits and credits, and those are split into a number of sub-series which share common traits (e.g. similar seasonality and/or trend). My question is, is it appropriate to   1. split the total series into component parts;   2. forecast each part in the most appropriate manner;   3. re-aggregate the parts into a complete projection? For example, I'd like to treat the following sub-series/parts differently: **Expenditures**   * Wages/Pension (annual percent increase, based upon contracts)   * Insurance (annual percent increase based upon prior increases)   * Supplies (annual percent increase based upon different prior increases) Then, add the individual sub-series back into a single forecast. Obviously, there would be similar but different sub-series/parts for revenue, which would eventually be combined for a total forecast.


title: Lower bound for non-parametric data, with a tolerance margin?

text: I have a set of data, which, from an initial test, is expected to be non- normal. My initial test shows a cluster of data gathered around 90, and another cluster gathered around 120. Neither cluster looks particularly normal either (but they both contain a very small amount of values, so it's hard to tell if they follow a precise distribution). From my knowledge of the problem, I'm expecting my real data to follow a similar pattern: a couple of clusters of data centered anywhere between 50 and 150. I'm trying to design an experiment to prove that the data will always be greater than zero. Intuitively, it should be reasonably easy to prove, given that all my data is generously above the limit. But I can't figure out a statistically sound way of doing this that doesn't require a prohibitive amount of samples. I've estimated, using Hahn & Meeker's formula as suggested in the NIST handbook in section 7.2.6.4, that in order to guarantee that 99.7% of my data is covered in the interval between the lowest and highest value, with 95% confidence, I'd need 1580 samples, which is not technically feasible for me. But that's not really what I'm trying to prove either. Is there a way I can take into account the large tolerance margin I have between my lowest sample, and the lower bound I'm trying to prove? Or am I inherently limited by the fact that I'm using to use non-parametric methods?


title: Account for group age differences

text: I have a dataset of two patient and one healthy control group which I would like to compare (using R) with respect to a continuous outcome variable (each subjects is measured once). However the groups differ in age. In the context of the dataset the age effect is not of interest. My questions are:   * What is the best way to account for the potentially confouding effect of age?    * From what I know one could do this by either adding it as a covariate in the linear model or matching the samples with respect to age via usage of the propensity score. Which method is recommended in which context?   * Is there an function to select subsamples of each group while matching groups with respect to age?


title: Is there any gold standard for modeling irregularly spaced time series?

text: In field of economics (I think) we have ARIMA and GARCH for regularly spaced time series and Poisson, Hawkes for modeling point processes, so how about attempts for modeling irregularly (unevenly) spaced time series - are there (at least) any common practices? (If you have some knowledge in this topic you can also expand the corresponding wiki article.) Edition (about missing values and irregular spaced time series) : Answer to @Lucas Reis comment. If gaps between measurements or realizations variable are spaced due to (for example) Poisson process there is no much room for this kind of regularization, but it exists simple procedure : `t(i)` is the i-th time index of variable x (i-th time of realization x), then define gaps between times of measurements as `g(i)=t(i)-t(i-1)`, then we discretize `g(i)` using constant `c`, `dg(i)=floor(g(i)/c` and create new time series with number of blank values between old observations from original time series `i` and `i+1` equal to dg(i), but the problem is that this procedure can easily produce time series with number of missing data much larger then number of observations, so the reasonable estimation of missing observations' values could be impossible and too large `c` delete "time structure/time dependence etc." of analysied problem (extreme case is given by taking `c>=max(floor(g(i)/c))` which simply collapse irregularly spaced time series into regularly spaced Edition2 (just for fun): Image accounting for missig values in irregularly spaced time series or even case of point process.


title: Average weighted to exclude a predictor variable?

text: I have over 200 datapoints in my set and I have related my predictors (mostly categorical) to a single continuous response variable in a multiple regression model. The model is to my satisfaction, and now I need help in analyzing one aspect of it. In my case, each datapoint represents an event, and my response variable is the number of attendees for that event. Increasing event frequency seems to decrease the number of attendees per event. However, so far increases in frequency have increased the OVERALL event attendance. So what I'd like to do is project the number of attendees for an entire week given various numbers for the total event count that week. Basically, I want to make a table with a column for the number of events that week, and another for number of total attendees for that week in order to select the optimal event frequency. I have a variable in my model for number of events held that week. I know that I can use its coefficient to determine the negative impact of holding one additional event in a given week, but that data isn't useful to me without a starting point to subtract it from. I'm looking for a way to predict the number of attendees for a given weekly frequency ignoring all other variables. Is there some kind of weighted average that will help me find this number?


title: How to define an informative prior in a varying intercepts (linear mixed) model

text: Gelman and Hill 2007 present a varying using WinBUGS and R: Simplifying the data description hugely, the data is from 85 counties, with log radon level measured in the basement and the first floor of houses in these counties (multiple houses in each county). y is log radon level, and x is floor (0 for basement, 1 for first floor, so treatment contrasts), county marks the id of the relevant county. This is data from Gelman and Hill's book. If this were pure R, we could do:               lmer(y~x+(1|county)).      Let's say the research question is: does being in the basement vs 1st floor affect radon emissions? I.e., is the slope b significantly different from 0 (it should be negative, because the higher up you are, the less the radon). I have WinBUGS code from Gelman and Hill's book (radon.1.bug), and modified code (radon.1a.bug), modified by me because I want an informative prior for the slope for county. I just used some rather large value for b as a prior, and a fairly low sd:               b ~ dnorm (10, .0001)      I first thought I could get away with making this one change; I expected the posterior estimate for b to be skewed towards my prior because the prior has such a low sd (reflects my strong belief that beta=10?). But this had absolutely no impact on the posterior! Can someone help me understand what I'm doing wrong? To run the model:   1. Download the following files from this site: http://www.stat.columbia.edu/~gelman/arm/examples/radon/ cty.dat, srrs2.dat, radon.1.bug (the WinBUGS model). Then download and run radon_setup.R to set up the data.   1. The code I used (Modified slightly from Gelman and Hill): radon.data <- list ("n", "J", "x", "y", "county") radon.inits <- function (){ list (a=rnorm(J), b=rnorm(1), mu.a=rnorm(1), sigma.y=runif(1), sigma.a=runif(1)) } radon.parameters <- c ("a", "b", "mu.a", "sigma.y", "sigma.a") library(R2WinBUGS) ## Specify the location of your WinBUGS directory: bugsdir<-"C:/Users/Shravan/Desktop/winbugs14/WinBUGS14" radon.1 <- bugs (radon.data, radon.inits, radon.parameters, "radon.1.bug", n.chains=3, n.iter=500,bugs.directory=bugsdir) plot (radon.1) radon.1.noburnin <- bugs (radon.data, radon.inits, radon.parameters, "radon.1.bug", n.chains=3, n.iter=500, n.burnin=0, bugs.directory=bugsdir) plot (radon.1) I modified radon.1.bug as follows (call this radon.1a.bug and replace in the above code to see what happens---it seems precisely nothing changes in the estimated slope):               ### Start WinBUGS model, call it radon.1a.bug:     model {       for (i in 1:n){         y[i] ~ dnorm (y.hat[i], tau.y)         y.hat[i] <- a[county[i]] + b*x[i]       }       #b ~ dnorm (0, .0001)  ## <- Change this to an "informative prior"       b ~  dnorm (10, .0001)  ## Shouldn't this influence the posterior?       tau.y <- pow(sigma.y, -2)       sigma.y ~ dunif (0, 100)            for (j in 1:J){         a[j] ~ dnorm (mu.a, tau.a)       }       mu.a ~ dnorm (0, .0001)       tau.a <- pow(sigma.a, -2)       sigma.a ~ dunif (0, 100)     }     ### End WinBUGS code.      I'd be grateful for any help, including pointing me to an example online or a book where an informative prior is used in the context of a normal-normal conjugate prior, in either a single t-test setting or a varying intercepts model of the sort described above.


title: Did I find a bug in the tseries or urca packages?

text: I'm trying the functions to check the cointegration of a matrix. I'm using **Phillips & Ouliaris Cointegration Test** The function in _tseries_ package is **po.test** and **ca.po** in _urca_ The results with **urca** are:               > ca.po(prices, demean='none')          ########################################      # Phillips and Ouliaris Unit Root Test #      ########################################           Test of type Pu      detrending of series none                Call:     lm(formula = z[, 1] ~ z[, -1] - 1)          Residuals:         Min      1Q  Median      3Q     Max      -7.4960 -0.2912  0.7116  1.4530  3.3962           Coefficients:             Estimate Std. Error t value Pr(>|t|)         z[, -1] 0.559705   0.004678   119.6   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 1.73 on 749 degrees of freedom     Multiple R-squared: 0.9503, Adjusted R-squared: 0.9502      F-statistic: 1.431e+04 on 1 and 749 DF,  p-value: < 2.2e-16                Value of test-statistic is: 12.9648           Critical values of Pu are:                       10pct    5pct    1pct     critical values 20.3933 25.9711 38.3413      The result with **tseries** are:               > po.test(prices, demean=FALSE)              Phillips-Ouliaris Cointegration Test          data:  prices      Phillips-Ouliaris standard = -25.6421, Truncation lag parameter = 7,     p-value = 0.01          Warning message:     In po.test(prices, demean = FALSE) : p-value smaller than printed p-value      As you can see I'm testing the same matrix (prices). How is it possible that **urca** tells there is **NO** cointegration and **tseries** : **YES**? Prices max it's a simple matrix with two columns (stock1 - stock2), take a look to an extract of that.               1  3.065448  5.244870     2  3.094924  5.806821     3  2.873858  5.647601     4  3.205457  6.190820     5  3.315990  6.453064     6  3.168612  6.865161     7  3.271777  7.230428      Thank you


title: What area of statistics is this example?

text: Let's say that there are a series of physical objects each with properties (round, red, soft, etc.) and each property affects the value of the object. Certain combinations can affect the value as well (round and red are valuable when found together, but not separate). Given a large list of objects and values, is there a way to ascertain the approximate contribution of value for each property and combination of properties? What is this called and where can I get more information? FWIW I am in the process of programming a module in Python. Also, I tagged this as homework. It isn't, but I couldn't find any other general tags. Thanks


title: Logistic regression with binary dependent and independent variables

text: Is it appropriate to do a logistic regression where both the dependent and independent variables are binary? for example the dependent variable is 0 and 1 and the predictors are contrast coded variables -1 and 1 ?


title: How to test whether percentage is significantly greater than 50% after averaging across three experiments?

text: ### Context: I ran 3 experiments comparing in each 2 techniques. In each experiment participants had simply to express a preference about which of the 2 techniques they thought was clearest to understand. 30 participants were involved, 10 in each experiment. No repeated measures were involved. Now I have a table of the results like so:                           Technique 1   Technique 2         EXP1:   80%           20%         EXP2:   90%           10%         EXP3:   60%           40%      I want to know if globally technique 1 is statistically significantly better than technique 2. ### Questions:   * Which analysis should I perform?    * How could this be done in R?


title: In SAS, is there a way to iterative refit a regression using different weights for observation with positive residual and negative residuals?

text: I have fitted a model using proc reg, say, using this               proc reg data = mydata;       model a = b;     run;      But in this particular application it is better to over-estimate than it is to under-estimate. So I actually want to refit this model this time using the residuals as weights so they an over-estimate is penalised more heavily. Is there a way to do that without writing my own macro?


title: observation and hypothesis testing

text: Can observation be used to collect data for hypothesis testing? Can we be confident with the power of the test in those cases? We know power of the test can be applied with hypothesis testing for data collected through sample.


title: Correcting for heteroskedasticity in Logit fixed-effects model

text: I found out that my panel data 'suffers' from heteroscedasticity by doing the test described on Stata FAQ, Testing for panel-level heteroskedasticity and autocorrelation. I understand that the $H_0$ for `lrtest` is homoskedasticity and I'm rejecting this. Please correct me if I'm wrong here. I'm carrying out a fixed effects logistic regression with `xtlogit` and I found that the `vce` option should be used to correct for heteroskedasticity. Still, I don't understand when to use `vce(bootstrap)` or `vce(jackknife)` ? What's the difference between the two? How can I decide which one to use? EDIT: Some info on the data size:                  Number of observations: 116304        Number of groups 2549        Obs per group: min = 1                       avg = 45.6273                       max = 20


title: Statistical analysis approach for seed germination and viability study

text: **Background:** I am a stakeholder in a research project involving seed germination and viability tests for a number of shrub species. The ultimate goals of this project are to determine best practices for future storage and use of these species' seed. To investigate this, seeds from a variety of accessions and ages, of numerous species, are subjected to various storage conditions (e.g.: air vs N2, high moisture vs. low moisture, three different storage temperatures, and different lengths of storage time [1-3 years]). Each experiment (e.g.: accession #1 of species A, collected in 2009, stored for 1 year with Air, and at 4 deg C, with low seed moisture), is tested for viability (%) and germination (%). Stratification measures have been applied to break seed dormancy before the experiment begins. **My problem is the following:** We have received preliminary data (2.5 years' worth of a 3 year study) for almost all species and accessions and tests, and the data are all over the place. A large part of this variability is likely due to the fact that shrub seed is not well-studied, and that the very deep dormancy is not well understood. The researchers have conducted ANOVAs on all of the experiments and their replicates, using initial germination and viability data as part of the comparison. However, this is where their analysis ends. They have identified many experiments that show statistical difference, however there are no further analyses to determine what factors may be causing this significance and having a strong impact on germination or viability. Because the data are so variable, I am desperately trying to determine a statistical approach to take to try to tease out which factor(s) may be showing the strongest influence on germination or viability. Each experiment combines a number of experimental factors (i.e. seed age, storage time, gas, temperature, moisture), so it would be great if I could understand if one of these factors should be investigated further. At first, I thought of performing a Principal Components analysis, because I envisioned some sort of ordination that would allow me to visualize how different groups emerge from the data, based on experimental factors. However, I later questioned this idea because all of the factors are controlled, and not measured outcomes of the experiment -- temperature wouldn't affect moisture level or seed age, etc., so I couldn't justify this analysis. **My question is this:** Is there some way to analyze these data? There are a handful of experimental variables that are controlled, but only one or two response variables (germination or viability over time). I thought that I could think of a strategy to gain more insights into the data, but I am completely lost.


title: Estimation parameters using Generalized Method of Moments and Generalized Empirical Likelihood with R

text: I have two density functions: $h_{1}(x)$, $h_{2}(x)$. I want to estimate parameters from one dataset using $ h_{1} $ and $ h_{2} $ andthe GMM estimator or the GEL : from the article "Computing Generalized Method of Moments and Generalized Empirical Likelihood with R" I want to estimate the parameters : $ \delta $, $\beta $ $ \mu $ and $ \sigma $. $h_{1}(x)=\dfrac{1}{\sigma \left( \frac{1}{6}-\frac{\delta}{2}+ \frac{\delta^{2}}{2} + \beta\right) }\left[ \frac{\left[ \Phi^{2}\left( \frac{x-\mu}{\sigma}\right)\right] }{2}-\delta \Phi\left( \frac{x-\mu}{\sigma}\right)+\left( \frac{\delta^{2}}{2} + \beta\right) \right] \phi\left( \frac{x-\mu}{\sigma}\right)$ Where $ \Phi $ and $ \phi $ are the cumulative and the density function of standard normal distribution. $h_{2}(x)=\dfrac{1}{P(1)} \left[ \dfrac{\left( e^{\frac{x-\mu}{s}}\right)^{2} \left( \frac{1}{2}-\delta+\frac{\delta^{2}}{2}+\beta\right) - \beta e^{\frac{x-\mu}{s}}+\left(\frac{\delta^{2}}{2} + \beta\right)}{s\left( 1+e^{\frac{x-\mu}{s}}\right)^{4}} \right].$ So my problem is how to do it with R, I am beginner in R but it is the only software that can deal with this method of estimation. I want to use GMM and GEL as in this article.


title: Inequality based conditional cumulative sum

text: I have been trying to do code this: For each $\left( x,z\right) \in\left( X,Z\right):$ $r(x,z)=\sum_{i=1}^{n}Y_{i}1\{X_{i}<=x\}1\left\{ Z_{i}<=z\right\}$ So far, the best way I came up to do this is by using a loop. Here an example:               y=rnorm(10)     x=c(1,1,1,2,2,2,3,3,3,4)     z=c(5,5,6,6,7,7,8,8,9,9)     data=data.frame(y,x,z)         s=rep(NA,length(unique(x))*length(unique(z)))     dim(s)=c(length(unique(x)),length(unique(z)))     for (i in 1:length(unique(x))){     for (j in 1:length(unique(z))){        s[i,j]=sum(y*as.numeric((x<=unique(x)[i]))*                     as.numeric((z<=unique(z)[j])))     }     }      The output is OK, but given that once the dimension of x and z grows, this becomes quite slow.Since, for a given z, this looks like a conditional cumulative sum, I am 100% sure that there is a more efficient way of doing this, without the loop. Would any of you have any suggestion? If I didn't have z, I know I could use data.table:                s=data[order(x)][,lapply(.SD, sum),by=c("x"), .SDcols=c("y")]       s=s[,lapply(.SD, cumsum), .SDcols=c("y")]      but with more than one index (x and z, not just x) I was not able to formulate the program.


title: Why Bayesian Networks? Why Markov Networks?

text: Undirected Graphical Models are usually known as "Markov Networks" and Directed Graphical Models are known as "Bayesian Networks". This is naming is not very clear to me. It might be because of some historical background, though I am just curious to know if there is any concrete reason for this naming.


title: Confidence intervals vs sample size?

text: I am totally new to stats and the field of confidence intervals. So this might be very trivial or even sound stupid. I would appreciate if you could help me understand or point me to some literature/text/blog that explains this better. I see on various news sites like CNN, Fox news, Politico etc about their polls regarding the US Presidential race 2012. Each agency conducts some polls and reports some statistics of the form: CNN: The popularity of Obama is X% with margin of error +/- x1%. Sample size 600. FOX: The popularity of Obama is Y% with margin of error +/- y1%. Sample size 800. XYZ: The popularity of Obama is Z% with margin of error +/- z1%. Sample size 300. Here are my doubts:   1. How do I decide which one to trust? Should it be based on the confidence interval, or should I assume that since Fox has a larger sample size, it's estimate is more reliable? Is there an implicit relationship between confidence itnervals and sample size such that specifying one obviates the need to specify the other?   2. Can I determine standard deviation from confidence intervals? If so, is it valid always or valid only for certain distributions (like Gaussian)?   3. Is there a way I can "merge" or "combine" the above three estimates and obtain my own estimate along with confidence intervals? What sample size should I claim in that case? I have mentioned CNN/Fox only to better explain my example. I have no intention to start a Democrats vs Republicans debate here. Please help me understand the issues that I have raised.


title: How do you find causal relationships in data?

text: Lets say I have a table with columns "A", "B" Is there a statistical method to determine if "A" causes "B" to happen? One can't really use Pearson's r, because:   * it only tests the correlation between values   * correlation is not causation    * Pearson's r can only correlate linear relationships So what other options do I have here?


title: Logistic regression assumptions for a model with many binary independent variables

text: I am working on developing a logistic regression model that uses qualitative variables only ($n=990$). My remit is to define the equation that can identify the most relevant characteristics of a survey respondent that is favorable towards Company X. The proposed equation is similar to the following: \begin{align} {\rm Fav} = &a*{\rm age} + b*{\rm CoAware} + c*{\rm IssueAware} + d*{\rm readnewspaper} + e*{\rm region} + \\ &f*{\rm income}\ldots \end{align} The dependent variable is "Company Favorability" (0 = Unfavorable/Neither | 1 = Favorable). There are currently 25 independent variables, 20 of which are binary IVs that range from highly correlated to the DV (awareness of Company) to not significant (gender). I also have 5 categorical variables that indicate region of the country, age (in categories), party affiliation, income level, and education. I am almost certain that I need to use a logistic regression model for this approach. However, when I test my assumptions, I have having a very difficult time proving a linear relationship between the dichotomous independent variables and the logit transformation of the DV. My other problem is that, I am somewhat overwhelmed by possible interaction effects. There are 34 possible options using 25 variables - leading me to over 50 million possible combinations. I have three questions:   1. Is there a better method to model with a binary dependent variable?    2. Am I missing something in the assumptions? (ie: Do I need to indeed prove the linear relationship if all of my variables are dichotomous)   3. Would it be better to approach this by looking at multicollinearity first, to reduce the number of variables overall, and then look at linear relationships with the logit of the DV?


title: Identifying when fixed effects have changed in a GLMM

text: I have a dataset spanning several years, that I have analysed year by year (using the same glmm model code) to determine whether the random effect (random intercept) has altered over time. I have found that it has altered. Now I am concerned that the fixed effects have also altered, which may mean that I can't compare the random effects year by year. To determine this, I have obtained confidence intervals for all the fixed effects for each year and compared the intervals. My thinking is that if the confidence intervals for a fixed effect for each year all overlap then that fixed effect can be assumed to be the same - is this faulty reasoning ? Is there a better way ? If this is a reasonable approach then does it imply that if the confidence intervals do not overlap then the fixed effect $has$ changed ?


title: Best practice when merging waves from non-panel survey

text: This is more of a data management question. When merging several waves of a survey where the same questions were asked in multiple years -- but each time of a new sample -- should the data from a question be merged into the same variable or separate variables? Since each wave was administered to a new sample, the unique (person) identifier would be new for each wave and there is no risk of analysing answers from multiple waves as though it were a panel dataset. But might using a single variable (say q1 rather than separate w1_q1 w2_q1 w3_q1) be confusing or helpful?


title: Iterative or Lazy Reservoir Sampling

text: I'm fairly well acquainted with using Reservoir Sampling to sample from a set of undetermined length in a single pass over the data. One limitation of this approach, in my mind, is that it still requires a pass over the entire data set before any results can be returned. Conceptually this makes sense, since one has to allow items in the entirety of the sequence the opportunity to replace previously encountered items to achieve a uniform sample. Is there a way to be able to yield some random results **before** the entire sequence has been evaluated? I'm thinking of the kind of lazy approach that would fit well with python's great itertools library. Perhaps this could be done within some given error tolerance? I'd appreciate any sort of feedback on this idea!